SEISMIC LIQUEFACTION ANALYSIS OF CAPITAL REGION OF ANDHRA PRADESH STATE , INDIA

Liquefaction is a phenomenon happens in a loose, fully saturated cohesionless soil in undrained condition subjected to cyclic loading. During liquefaction of the soil lost its shear strength when the mean effective stress is made equal to zero due to the progressively increasing excess pore water pressure. Liquefaction may cause failure of foundations, resulting in collapse of structure, even if the structure is designed as an earthquake-resistant. Liquefaction depends on characteristics of subsurface soil. Amaravathi is a new capital of Andhra Pradesh State, India. The construction activities in the capital region are swiftly increasing. It is essential that the new structures constructing in capital should be assessed for liquefaction susceptibility. In the present investigation an attempt has been made to assess the liquefaction susceptibility of various sites in the capital of Andhra Pradesh State, India. The liquefaction analysis is carried out by using simplified method which mainly relies on Standard Penetration Test (SPT) value.


INTRODUCTION
Earthquakes are most powerful natural disasters which are unavoidable. The hazards associated to earthquakes are referred to as seismic hazards. During an earthquake there is release of energy which reaches to the ground surface and to the structures by means of seismic waves. One of the major causes of destruction during an earthquake is the loss of strength & stiffness of cohesionless soils. This phenomenon called liquefaction occurs mainly in loose & saturated sand. When an earthquake shakes loose saturated sand, the grain structure of soil tends to consolidate into more compact packing. The soil liquefaction depends on the magnitude of earthquake, intensity & duration of ground motion, the distance from the source of the earthquake, site specific conditions, ground acceleration, type of soil and thickness of the soil deposit, relative density, grain size distribution, fines content, plasticity of fines, degree of saturation, confining pressure, permeability characteristics of soil layer, position & fluctuations of the ground water table [1,2]. The purpose of the study is to evaluate the liquefaction susceptibility of various locations in the capital of Andhra Pradesh state, India using penetration resistance value from standard penetration test after necessary corrections. Firstly Cyclic Shear Stress Ratio (CSR) that would be induced due to earthquake was computed. In calculating CSR, the peak horizontal ground = . ( Where, amax is the Peak ground acceleration g is the acceleration due to gravity = z is the Zone factor Seismic zoning map of India prepared based on the peak ground acceleration (PGA) induced by the maximum considered earthquake [3].
In the above two equations Z is the depth of the soil stratum The CYCLIC SHEAR RESISTANCE RATIO (CRR) is calculated from the following equation: Where (N1)60CS is the corrected SPT value including correction for fines

Factor of Safety against Liquefaction (FOS) is ratio of CYCLIC SHEAR RESISTANCE RATIO
(CRR) to CYCLIC SHEAR STRESS RATIO (CSR). 2.56 (7) In which M is the Magnitude of the Earthquake If the value of Factor of Safety against Liquefaction is less than or equal to 1, the soil is susceptible to liquefaction (L). If the value of Factor of Safety against Liquefaction is greater than 1, the soil is not susceptible to liquefaction (NL). Standard Penetration Test (SPT) was conducted as per the guidelines of [13]. The SPT is carried out in drilled boreholes, by driving a standard 'split spoon' sampler using repeated blows with a 63.5kg hammer falling through 750mm. The bore holes have been drilled using rotary hydraulic drilling of 150mm diameter up to the rock depth. The hammer is dropped on the rod head at the top of the borehole, and the rod head is connected to the split spoon by rods. The split spoon is lowered to the bottom of the hole, and is then driven for a depth of 450mm, and the blows are counted normally for each 150mm of penetration. The penetration resistance (N) is the number of blows required to drive the split spoon for the last 300mm of penetration. The penetration resistance during the first 150 mm of penetration is ignored. The 'N' values measured in the field using SPT procedure have been corrected for various corrections recommenced for evaluating the seismic borehole characteristics of soil.
First, corrected 'N' value i.e., (N1)60 are obtained using the following equation: ( 1 ) 60 = Then corrected 'N' values (N) is further corrected for fines content based on the revised boundary curves derived by [14] as described below: The N value for soil shall be corrected for overburden is extracted from [11].
Where is the effective overburden pressure.
Correction for Hammer Effect [CER] can be taken as follows: For Doughnut hammer : 0.5 to 1.0 Liquefaction, in the past, was primarily associated with medium to fine grained saturated cohesion less soils and soils with fines were considered non-liquefiable. [13] studied the liquefaction behaviour of silts and silt clay mixers over a range of plasticity values of interest by conducting cyclic tri axial tests on undisturbed as well as reconstituted samples and their behaviour was compared with that of sand. Saturated silts with plastic fines were found to behave differently from sands both with respect to rate of development of pore water pressure and axial deformations. Later on it was found by several investigators that certain soils with fines may be susceptible to liquefaction. Δ(N 1 ) 60 = e Where is the fines content Corrected SPT value including correction for fines [(N1)60CS] is given by (N 1 ) 60CS = (N 1 ) 60 + Δ(N 1 ) 60 (11) Tab. 1

-Liquefaction Analysis of Mandam Site for an Earthquake Magnitude 4
Liquefaction analysis has been conducted for other earth quake magnitudes also. Analysis carried in various locations of capital region of Andhra Pradesh State, India. Depth versus factor of safety against liquefaction for different Magnitudes of Earthquake at various locations of capital of Andhra Pradesh State is shown in Figures 3 to 10. It was observed from Figures 4 to 11, the Factor of Safety against Liquefaction is greater than 1 for Earthquake Magnitudes of 4 and 5 irrespective of depth for various locations. It was also noticed that when Magnitude of earthquake 6 and above, many of sites in the capital region of Andhra Pradesh State, India was prone to Liquefaction. At some locations (as depicted in Figure 7, 8 and 10) higher factor of safety against liquefaction was obtained due to presence of hard gravel (SPT value is more than 50) at lower depths. Depth   The number of sites suceptible to liquefaction with increase in earthquake maginitude is plotted in the Figure 11. By observing the Figure 11, it was noted that the majority of the locations in the region are liquefiable when the magnitude of quake exceeds 5. It was unequivocally distinguished that no site has indicated liquefaction vulnerability for a seismic quake for maginitude of 4. It can be inferred that the greater part of the capital region of Andhra Pradesh may not liquefiable when a light seismic quake happens.

CONCLUSION
Liquefaction susceptibility assessment based on empirical approach with SPT value (N value) was carried out at various places in capital of Andhra Pradesh state, India. The following conclusions were drawn from the analysis.
1. Liquefaction analysis should be carried for the all the structures and mainly for structures of national importance. 2. Many of the liquefaction susceptibility studies need to be considered in early stages of planning and design in order to select the most appropriate sites and also to improve sites for mitigating liquefaction susceptibility. 3. Liquefaction susceptibility assessment based SPT value is quite feasible. Since SPT was most widely used method soil exploration in India. 4. Most of the sites in the capital region of Andhra Pradesh are susceptible to liquefaction at the magnitude of earthquake 6 and above. 5. Majority of the sites of capital of Andhra Pradesh may not susceptible to liquefaction when a light earthquake happens. If a strong earthquake happens almost all the areas in the capital of Andhra Pradesh are susceptible to liquefaction and there is a chance of huge loss for life and property. 6. The findings would help the designers in taking suitable decisions for design of foundations, resistant to liquefaction and to adopt appropriate ground improvement techniques for rapidly developing capital of Andhra Pradesh. [1] Kramer, S.L. "Geotechnical Earthquake Engineering", Pearson Education (Singapore) Private Ltd., New Delhi, 2003.

INTRODUCTION
Bridge deflection is an important basis to evaluate bridge health as it directly reveals the rigidity, stability, bearing capability and earthquake resistance of the bridge [1]. Bridge dynamic deflection can reflect the load impact coefficient and the internal force distribution of the structure [2]. It is therefore important to monitor the dynamic deflection of a bridge with a traffic light as it confronts heavier vehicle dynamic loads than others [3].
At present, there are many methods to monitor bridge deflection such as the precise levelling, the total station, the measurement robot, the laser speckle, the photographic imaging, the inclinometer method, GPS (Global Positioning System), the inertial measurement, the tension line method, the pipe deflection monitoring method and the photogrammetric techniques. The precise levelling can monitor the static bridge deflection with high precision, it is however out of monitoring the dynamic bridge deflection [4]. The total station can be used to monitor the deflection of a longspan bridge. It is however challenging to monitor the bridge dynamic deflection [5]. Although the measurement robot can monitor the instantaneous coordinates of a single point, it only can monitor the short periodic deflection of the bridge as it continuously monitors the same point two times with a periodic [6]. The laser speckle and the photographic imaging can monitor the dynamic deflection of a point on the bridge without monitoring the dynamic global deflection of a bridge [7,8]. The inclinometer can monitor the bridge dynamic deflection with high precision, but it requires the installation shaft parallel to the bridge axis which is difficult to carry out in the field [9]. GPS can be used to monitor a long-span bridge as it has a low accuracy in monitoring the dynamic deflection of a bridge [10]. The inertial measurement has a high-resolution ability, but it is ineffective in low frequency [11]. The tension line method cannot monitor the global dynamic deflection of a bridge as one tension line only can monitor the deflection of one point [12]. The pipe method can monitor the static global deflection of the bridge, but it is not flexible in monitoring the dynamic deflection of a bridge yet [13]. Photogrammetric techniques [14,15] can monitor the global deflection of a bridge, and it has advantage in non-contact measuring. But it cannot monitor the dynamic global deflection of a bridge as it uses two cameras or takes images from at least two different positions by one camera [16].
As such, a specific method is required to monitor the dynamic global deflection of a bridge and warn the possible danger of the bridge in time. MDP offers the potential to solve this problem [17,18]. MDP, combining close range photogrammetric technique [19][20][21] with information technology [22,23], can monitor dynamic global deflection of a bridge and obtain image sequences of a bridge as it continuously monitors a bridge by a single digital camera. Although it has not been as popular in bridge structures as the photogrammetric techniques, many pioneering applications in this field have proved its increasing capability [24,25].
It is clear that in most of the previous studies MDP had been used in monitoring bridge deflection. However, they did not consider problems such as digital camera parallaxes caused by the environment and the photographing direction being un-perpendicular to the bridge plan.
The aim of this study is to propose the PST-TBP (photographing scale transformation-time baseline parallax) method in MDP to monitor the dynamic global deflection of a bridge to grasp the deflection characteristics of the bridge caused by vehicle dynamic load and MDPS (monocular digital photography system) is used to depict the deflection trend curves of the bridge in real time to assess bridge health on site and warn the possible danger. − 1 + 2 + 3 + 4 9 + 10 + 11 +1 = 0 − 5 + 6 + 7 + 8 9 + 10 + 11 +1 = 0 where x and z are image plane coordinates of deformation points without errors, X, Y and Z are the space coordinates of the correspondence deformation points, and Li (i =1, 2,3,…,11) are the functions of the exterior and interior parameters of a digital camera.
The error (2) can be obtained by linearizing Equation (1): where P1 and P2 are the weight matrices of image point observations and the reference point observations, respectively, ∆ is the correction matrix of the reference points, N0 is the coefficient matrix of ∆ ，M is the coefficient matrix of ∆ , and 1 is the constant matrix of image point observations.
In Table 1, the calculation distances of U0-U2, U1-U3, and U2-U4 were obtained by the DLT method. Measurement distances of U0-U2, U1-U3, and U2-U4 were seen as precise values. The measurement errors were obtained by differencing the calculation distance with the corresponding measurement distance. The maximal measurement error of the digital camera is within 1mm which meets the accuracy requirements of deformation monitoring [27].

A principle of photographic scale transformation
The photographic scale of somewhere always changes along the photographing distance (from the position to the photography centre) [28,29]. Figure 1 shows a schematic diagram of a CCD (Charge Coupled Device) camera capturing images at different photographing distances H3 and H4. H1 is the focal length of a CCD camera, H2 is the distance between the optical origin (o) and the front end of CCD camera, D1 on reference plane and D2 on object plane are the real-world length formed by the view field of the CCD camera at photographing distances H3 and H4 respectively, and N is the maximal pixel number in a horizontal scan line of an image plane, which is fixed and known as a priori irrelevant of the photographing distances.   Figure 1, the relationship between pixel counts and distances can be described by: In general, H3 and H4 are meter-sized, while H2 is centimetre-sized. Assume that H2 can be ignored when the camera is far from the bridge, Equation (3) can be expressed as: From Equation (4), we have: Assume that M1 and M2 are the photographing scale of the reference plane and the object plane, respectively. According to Equation (5), we have: Namely, where ∆ is the photographing scale transformation coefficient, and ∆ = 4 3

Photographing scale transformation-time baseline parallax method
Assume that there are no errors in the measurement, the horizontal and vertical displacements of deformation point on object plane based on the time baseline method are given by: where and are the horizontal and vertical deformation of a deformation point on the object plane, ∆ and ∆ are the horizontal and vertical displacements of the corresponding deformation point on the image plane, M is the photographic scale on the reference plane, and is the photographing scale on the object plane. Note that ∆ and ∆ are with parallax errors.
However, a time baseline parallax method requires the photographing direction perpendicular to the bridge plane and the camera constantly when a single-digital camera is used to monitor the bridge. It is difficult to carry out [30]. This study proposes the PST-TBP (photographing scale transformation-time baseline parallax) method to solve these problems.
The PST-TBP method consists of three steps. Firstly, the TBP method is used to get the displacements on the reference plane of deformation points. These displacements are with the parallax errors caused by the change of intrinsic and extrinsic parameters of a digital camera. Secondly, reference points are used to match a zero image with the successive images to eliminate the parallax errors. And the corrected displacements on the reference plane are obtained. Lastly, the real displacements on the object plane of deformation points are equal to the corrected displacements on the reference plane multiplied by the photographing scale transformation coefficient. The details are as follows: The first step, when a point on an object plane moves from A to B (Figure 2), its horizontal and vertical displacements on reference plane are given by: where ∆ and ∆ are the horizontal and vertical displacements on reference plane of a deformation point, ∆ and ∆ are the horizontal and vertical displacements on image plane of deformation point, and M is the photographing scale on the reference plane. Note that ∆ and ∆ are with parallax errors.
In the second step, some points are laid at a stable position around the camera to form a reference plane which is perpendicular to the photographing direction, and the parallax is therefore eliminated through differencing a zero image with the successive images based on reference plane, respectively.
After correcting the displacements on the image plane of a deformation point based on the parallaxes of the reference plane, the corrected displacements on the image plane of a deformation point is obtained: where (∆ ′ , ∆ ′ ) and (∆ , ∆ ) are the corrected displacements and the measured displacements on the image plane of a deformation point, respectively.
Thus, we obtained the corrected deformation values on the reference plane of a deformation point: where(∆ ′, ∆ ′) are the corrected deformation values on the reference plane of a deformation point.

BRIDGE TEST
Before the test, we set the total station and three digital cameras on the specified position of the south side of the Xiaoqing river. To facilitate the stress analysis of this bridge and the field condition, six deformation point targets labelled as U0-U5 were set on the bridge's superstructure, and seven deformation point targets labelled as U6-U12 were evenly set on the bridge's major deck. Reference points labelled as C0-C1, forming the reference system, were set near the digital camera whose view is in Figure 4. The reference system, used to match a zero image with the successive image, is perpendicular to the photographing direction. July 16, 2009) In the bridge test, the SONY350 cameras were held constantly as much as possible. The cameras used in this test can capture the instantaneous deflection of a bridge in 1‰ second and shoot this bridge seven times in one second. The test process was described as follows:

Fig. 4 -Field test on the Fenghuangshan road bridge (taken on
(1) In the test, one excavator moved on the deck at a speed of 20km/h from the north to the south. The cameras were used to shoot the bridge to produce the zero images before the excavator on the bridge.
(2) The cameras were used to shoot the bridge every two seconds when the excavator was moving on the bridge. Finally, ten successive images were produced.
(3) The total station was used to obtain the spatial coordinates of reference points and deformation points after the test. Table 2 shows the two-dimension coordinates of selected monitoring points and camera as the elevation can be ignored in the study.

Photographing scales of deformation points
Based on the photographing scale transformation, the photographing scales of U6, U9, and U12 can be expressed as: where 6 ， 9 and 12 are the photographing scales of U6, U9 and U12 after the photographing scale transformation，respectively; 0 is the photographing scale on the reference plane; ∆ 6 ，∆ 9 ∆ 12 are the photographing scale transformation coefficients of U6, U9, and U12, relative to the reference plane. 12 can be expressed as: where OA, OB, OC, and OD are the photographing distances of the reference plane, U12, U9, and U6, relative to the camera. They are detailed in Figure 5. The photographing direction is perpendicular to the reference system, Line4, Line5, and Line6. Point A is the projection of Line 2 in reference system. Points B, C, and D are the intersections of Line 4 and Line 2, Line 5 and Line 2, Line 6 and Line 2, respectively. OB is the projection of Line1 in Line2 and OD is the projection of Line 3 in Line2.
Similarly, the approximate photographing scales of C0, C1, U7, U8, U10, and U11 were obtained. Table 3 shows the photographing scales of these reference points and deformation points on the bridge deck.

Measurement accuracy
In theory, reference points did not move during the test and their displacements were zero. However, the displacements of these reference points obtained by the PST-TBP method were not zero. These displacement values of these reference points can therefore be used to represent the measurement accuracy. Table 4 shows that the maximal error is 1.65mm in monitoring the reference points.
As image matching is the key to this monocular digital photography, we measured the image coordinates of these deformation points on the zero image times to assess its image matching accuracy. Table 5 shows that the maximal error was 1 pixel, and the minimal error was 0 pixels. The average errors of U6-U12 were 0.1 pixels, 0.1 pixels, 0 pixels, 0.3 pixels, 0 pixels, 0.3 pixels and 0.2 pixels, respectively. This method reached a sub-pixel Image matching accuracy.

Analysis of bridge deflection trends
We calculated the pixel displacements (Table 6) and space displacements of some deformation points ( Table 7). The positive and negative values (Tables 6 and 7) represent the deformation point moving up and down, respectively. The deformation points on the bridge deck (U6-U12) were chosen as the aim of this paper is to study the bridge deflection trends caused by vehicle dynamic load. Table 7 shows that in Test 6 deformation point U7 developed the maximal deflection-55.34 mm which was within the allowed value-75mm (the allowed value = the bridge span/800, where L = 60 m in this study). Deflection curves (Figures 6 and 7) were also depicted in order to visually analyse the bridge deflection law caused by vehicle dynamic load. Figure 6 shows that the deflection of every position on the bridge's major structure was inelastic, all of which was within the bridge allowed value. Figure 7 shows that vehicle dynamic load results in the bridge moving down in spite of U7 moving up in test 5 and 7. In addition, the bridge global deflection curves fluctuated up and down like some sinusoidal-cosinusoidal curves. The bridge deflection was a parabola when the excavator moved to the bridge centre (Test 6). Especially, the deflection of every position on the bridge almost reached its maximal when the excavator moved to the bridge centre, which conforms to the deformation characteristics of the bridge caused by the external load. Note that it is impossible to monitor the bridge deflection with high accuracy when the camera is far from the bridge. But it can provide data support to study the deflection deformation law of the bridge caused by vehicle dynamic load. Although we cannot get the precise deflection of the bridge, the deflection trend of the bridge is also an important indicator to assess the bridge health. In the future, MDP & MR (monocular digital photography and measurement robot) would be a useful method to monitor the global deflection of the bridge as the measurement robot can monitor precise short periodic deflection of the bridge, and MDP can monitor the global deflection trend of the bridge.

CONCLUSION
This study used MDP based on the PST-TBP method to monitor the dynamic global deflection of a bridge. Before the bridge test, we first set SONY350 cameras in proper place and levelled them. Then, the reference system, formed by reference points, perpendicular to the photographing direction, was set near the selected camera. We produced a zero image using the cameras to shoot the bridge before the test and produced image sequences by shooting the bridge every two seconds when the excavator was moving on the bridge. Through the results, the following conclusions are obtained: (1) The PST-TBP method reaches a sub-pixel image matching accuracy. The average image matching error was within 0.3 pixels. And the maximal error was 1.65mm and 5.91 mm in monitoring the reference points and dynamic global deflection of a bridge, respectively. (2) Every position of the bridge almost reaches the maximal deflection when the excavator moves to the bridge centre. The maximal deflection of the bridge was 55.34mm which was within the bridge allowed value-75mm. (3) The deflection curves of the bridge fluctuated up and down like some sinusoidalcosinusoidal curves. The bridge deflection was a parabola when the excavator moved to the bridge centre. (4) Deflections of every position on the bridge were inelastic, all of which was within the bridge allowed. This indicates that the bridge is in good health.
The MDPS (monocular digital photography system) used in this study can monitor the dynamic global deflection of a bridge. Deflection trends of the bridge depicted by this system in real time show the bridge deflection visually, which is effective in warning the possible danger of the bridge. Although it cannot get the precise bridge deflection when the camera is far from the bridge, the dynamic deflection trend of the bridge is also an important indicator to assess the bridge health. In the future, MDP & MR (monocular digital photography and measurement robot) would be a useful method to monitor the bridge deflection as it can get precise short periodic deflection deformation of the bridge and the dynamic global deflection trend of the bridge. [11] Y. Xie, "Inertial Measurement Method for Railway Bridge Dynamic Deflection," Chinese Journal Ofentific Instrument, 1999. [12] J. F. Stanton, M. O. Eberhard, and P. J. Barr, "A weighted-stretched-wire system for monitoring deflections," Engineering Structures, vol. 25, pp. 347-357, 2003. [13] Y. Liu, Y. Deng, and C. S. Cai, "Deflection monitoring and assessment for a suspension bridge using a connected pipe system: a case study in China

INTRODUCTION
Water erosion is the most widespread form of soil degradation. Reducing the erosion is one of the major challenges in landscape management. Sediment transport from arable land into surface waters (streams, rivers, reservoirs) induces many significant problems in water bodies [1]. This problem was highly accelerated when the landscape patterns were destructed during technical improvement of agro technology and land consolidations. In Eastern Europe this effect was accelerated by collectivization of agriculture [2]. Based on research it is known that sediment transport can be reduced by inhibition of the surface runoff and proper agricultural management.
Rainfall simulation experiments are widely used as a method to study various flow and transport processes induced by rainfall [3]. They have been used on different land uses, slopes, scales, soils and climate conditions [4], [5], [6]. But most of the experiments have been driven by rather narrow objectives and therefore do not provide data for complex description of rainfall runoff and sediment transport processes. That is one of the actual scientific objectives nowadays -to describe the driving mechanisms of the physical processes in details and consequently to implement the plot scale data in optimization of the catchment scale parameters.
Surface runoff rate and sediment yield in the form of suspended solids concentration are standard variables observed in experiments oriented on soil erosion research with use of rainfall simulators. The review of small portable simulators used across Europe provides [7]. Due to the nature of the device the small simulators with watered area around 1 m 2 are more frequent. Larger simulators can be found as well [8], [9], [10], [11]. They are essential for a consequent mathematical modelling of both surface and subsurface runoff and related processes.
Together with eroded soil particles various nutrients are transported in the surface runoff. Phosphorus is the most significant and monitored particulate nutrient in the environment [12], [13]. It enters the water courses and is retained in the ponds and water dams where it impairs the environmental balance and water quality [14].
Another runoff component is the quick subsurface flow also known as hypodermic flow or interflow [15]. On the arable land it mostly occurs on the interface between the shallow ploughed topsoil and compacted subsoil. Since these two flow components may interact and lead to implications for both erosion and hydrologic response to the rainfall event, it is important to observe and model the subsurface processes as well.
Soil erosion intensity (in particular splash type) is decreased with the vegetation growth [16]. The soil surface is protected largely due to the plant's leaves, which absorb part of the kinetic energy of the raindrops. Farming operations also directly affect soil movement through activities such as tillage [17], root crop harvesting, and the trampling of soil and removal of vegetation by livestock. Combination of the current state of vegetation and farming operations influenced state of soil aggregates. Soil aggregates are less exposed to the direct impact of the raindrops and associated crumbling, which is considered as a trigger of the erosion process [18]. In order to evaluate the effect of the vegetation on particular soil erosion event it is therefore necessary to document the actual leaves extent. To this end two most common parameters are used: Canopy Cover and Leaf Area Index (LAI). Repetitions of the experiments in the same site during the growing season as well as comparing results from vegetated and bare experimental plots enable the vegetation cover and state of the soil surface effects to be assessed.
Together with detaching soil particles from the aggregates by the impact of raindrops also the surface runoff plays an important role in soil erosion, in particular its volume and flow velocity. Higher velocity induces larger shear stress and dragging forces to carried soil particles. From the known particle-size distribution of the eroded material basic physical principles of soil erosion can be verified [19].

Device and setup description
Mobile rainfall simulator operated by CTU in Prague was constructed in 2012 in cooperation with Research Institute for Soil and Water Conservation (VÚMOP). Its detailed description can be found in [20] and only basic information will be given here. The essential part of simulator is a folding boom with four nozzles by Spraying Systems -FullJet 40WSQ [10] controlled by electromagnetic valves. The boom consists of a dural steel framed structure and eight supporting telescopic legs. The boom can be folded out directly from the trailer and heaved by a pulley into desired height. Once the legs are joined it is completely standalone and can be detached and moved independently, constrained only by the length of the control wires and water supplying hosepipes. The whole construction was designed to be easily and quickly set out and packed back again without any interfering with the plot marked for the simulation. In order to prevent the wind from affecting the simulation the supporting construction is covered with the tarpaulin.
The device includes a 1 m 3 tank, portable generator, pump producing a constant pressure with minimum delivery 40 l/min and a control unit for managing the rainfall intensity by triggering the electromagnetic valves. The electromagnetic valves trigger the separate pairs of nozzles. Their opening and closing in an arbitrary time interval produces the desired intensity of simulated rainfall. Coupling the nozzles into pairs reduces hydraulic shocks in the water distribution system and therefore helps maintaining a uniform spatial distribution of water over the experimental plot.
Four nozzles FullJet 40WSQ are positioned 2.65 m above the ground and 2.4 m apart. Type of the nozzles as well as the position were selected after many calibration tests aimed at reaching as uniform spatial distribution and droplet characteristics as similar to natural raindrops as possible. These tests were carried out on the plot 2.1 x 2.8 m centred under the first nozzle (for exact position see Figure 1), assuming that the rest of the area is symmetric. Final configuration yielded the Christiansen's uniformity index [21] of 80% and more for all rainfall scenarios commonly used in the field simulations.
This configuration and water pressure 0.7 bar leads to the total watered area 2.4 x 9.6 m. It can be used for variable setup of experimental plots, although limited by decreasing rainfall intensity towards the edges of the watered area. In CTU simulations there are two separate plots used simultaneously, one with dimensions 1 x 1 m and another 2 x 8 m with the longer edge along the slope, as illustrated in Figure 1.

Rainfall characteristics
The most crucial variable influencing the course of the simulations are the rainfall intensity (spatial distribution and average over the plot) and droplet characteristics. Droplet size and impact velocity were analysed by a disdrometer, spatial distribution of intensity by a network of small buckets [22]. Average drop size (d50) reaches 1.75 mm for 0.7 bar pressure. Christiansen uniformity index CU [23] reaches 85%.
In order to maintain uniform spatial distribution of water over the plots, the experimental setup was designed with rather large overspray of the experimental plots. Therefore only part of the sprinkled water falls onto the plots and a final intensity needs to be checked in each simulation.
First the outflow rate at the nozzles is checked by collecting and weighing the sprinkled water in a four minutes test, as it differs slightly depending on the exact water pressure in the system. Next a reference measurement is performed at the beginning of every campaign on the fallow plots covered with an impervious tarpaulin. The rainfall intensity is derived from the total discharge after the steady state is reached. From the acquired values from a whole set of campaigns a calibration curve describing the relationship between the outflow at the nozzles and resulting intensity onto the plots is being derived and consecutively refined. The curve is then used for rainfall intensity determination in simulations on the vegetated plots.
The precipitation intensity is checked in 10 a step during every simulation run with the tipping bucket rain gauge with collecting area 200 mm 2 and resolution of 0.2 mm set on a fixed location above the upper end of the larger experimental plot ( Figure 1). Its output represents point values which cannot be used for deriving mean rainfall intensity over the plot but is valuable for verification of the temporal uniformity. Second continuously recorded variable is the water pressure at the first nozzle, which is monitored with the digital pressure sensor. Recording interval of 4 s is used due to the short intervals of valves operation. Acquired data serves for detecting eventual pressure changes in the water distribution system, which could affect the rainfall intensity and bias the simulation outcomes. As illustrated in Figure 2 lapses in rainfall intensity or even stoppage of the simulator occur and have a fast response in observed runoff. Data from both rain gauge and pressure sensor is necessary for correct water budget balancing and explanation of induced variability in observed variables. Moreover it is essential in experiments with variable rainfall intensity when reproducing natural rainfall events.

Surface runoff and plot scale effect
Time to initiation of the surface runoff is determined visually by observing the discharge from the plot in the lower outlet. The runoff rate is measured in situ gravimetrically for variable time intervals. Within first 10 minutes after the runoff initiation the interval is 1 minute, another 10 minutes the runoff rate measurements are performed in 2 minute intervals and for the rest of the simulation (until 1 hour after the runoff initiation or reaching steady state) the interval is increased to 5 minutes. The purpose of this arrangement is to capture the dynamics of the runoff rate and suspended solids concentration. They change considerably within the initial phase and tend to steady in the course of the simulation, as shown in Figure 2.
The runoff needs to be measured independently for each of the two parallel plots. However the discharges are measured in the same intervals and samples are taken at the same relative time with respect to corresponding times of initiation. This way the runoff dynamics in both plots (1m 2 and 16 m2) can be compared.

Soil water regime and shallow subsurface flow
The runoff monitoring is accompanied with soil water regime monitoring. The soil water content is measured in two ways. Undisturbed 100 cm 3 soil samples are taken both before and after the simulation and later they are analysed in the lab. The soil water regime in the top 5 cm of the soil profile is continuously monitored with the ThetaProbe ML2x and recorded in 1 minute interval. Additional probes are installed in the depths of 30 cm and 50 cm below the surface.
To monitor shallow subsurface runoff a trench is being excavated along the lower edge of the larger experimental plot. It is 2 m long and 0.6 m deep. At the bottom of the trench a drainage pipe with diameter of 100 mm is installed. Both the bottom and the walls of the trench from the topsoil / subsoil divide below are covered with an impermeable plastic foil so that only shallow runoff from the topsoil reaches the drainage pipe. The trench is partly filled with the gravel and enclosed by the topsoil. The drainage pipe is connected through the collection pipes to a tipping bucket, where the subsurface runoff is continuously measured or sampled. The collection system was designed in a way which does not require removing prior to each agrotechnical operation.

Maximum flow velocity
For the flow velocity determination, two methods are used in parallel. For the maximum flow velocity which corresponds to the rillflow in the preferential flow paths the progress of the dye tracer Brilliant Blue FCF is observed. This method is especially used on the cultivated fallow plots. From the top of the plot the time intervals are measured for the tracer to travel each 1 m distance. In order to highlight the head of the moving tracer an extra dose of the dye is added at the beginning of each 1 m section. This measurement is repeated three times for each simulation after reaching the steady outflow and ideally close to the moments when the suspended soils samples are taken. The timing is important for the verification of the relationships between flow velocity and the grain size of suspended particles.

Average flow velocity
The dye tracer method is partially subjective due to the need of visual tracing and it is hard to use in case the surface is covered with full-grown vegetation. Besides the maximum flow velocity the mean velocity represents the process of sediment transport. This variable is obtained by measuring the electric conductivity after the application of potassium bromide (KBr) or sodium chloride (NaCl) solution. It is usually applied subsequently in 2, 4 and 6 meters from the lower end of the plot, so it is possible to derive the velocities in separate sections of the sloping plot. Using a conductometer (WTW Cond 3310) the beginning of the rise and the time of the peak of the conductivity wave is observed. Due to normally very long recession limb of the conductivity wave this measurement is being stopped once a distinct lapse in conductivity is observed. Beside the difficult manipulation (conflict with other measurements) another disadvantage of this method is a slow yet unwanted soil salinization which might in the long term lead to the inhibition of vegetation growth and impairing the soil fauna in the location.

Suspended solids
In order to evaluate the dynamics of suspended solids concentration the runoff samples are taken in the same time intervals as in the runoff rate measurements. The samples are weighed in the field and later analysed in the laboratory according to the standards [24]. Suspended solids are filtered from the solution, dried and their weight is converted into concentrations. Final sediment yield results from relation of concentrations to the runoff rates in respective time intervals of sampling.
Along with the samples for concentration analysis another set of samples is taken for determining the particle-size distribution. Since such determination is rather time-consuming, these samples are taken in 10 minute interval only. Based on results of particle-size analysis the changes in sediment composition during the erosion event may be monitored. When related to the runoff velocity the particle-size may also be used for determination of friction and dragging forces in the overland flow.

Phosphorus transport
Concerning the chemicals the phosphorus is monitored. Phosphorus concentration analyses are carried out on mixed samples composed of several 50 ml of runoff taken at different times. Resulting 5 samples represent three 5 minute and two 15 minute intervals of runoff. There are two basic forms of phosphorus in surface runoff: particle-bound and dissolved phosphorus. Mixed samples are stabilized in the field with highly concentrated HNO2 (0.5 ml HNO2 in 100 ml of sampled water) in order to avoid dissolved phosphorus decay. Later in the laboratory two types of samples are analysed using ICP-OES, Inductively coupled plasma optical emission spectrometry [25]. In the original unfiltered samples the total phosphorus concentration is determined. For dissolved phosphorus determination the suspended solids are removed from the samples using 0.45 µm filter. The particle-bound phosphorus is determined then as the dissolved phosphorus subtracted from the total.

Leaf Area Index
LAI is measured within the LAI-2000 Plant Canopy Analyzer. The device is based on comparison of the light intensities above the vegetation and near the soil surface. Advantage of this method is that it is fast and enables repetitive or multiple measurements for obtaining a representative value. The measurement is always performed on the vegetated experimental plot and in the open crop field next to it. The value from the latter site serves as a reference value for detecting significant variance in the vegetation growth in the experimental plot due to previous simulations.

Canopy Cover
The measurement of CC is performed by sensing of the surface with the camera and analysing the images using a supervised classification in GIS software. The analysed plot is delimited with a metal squared frame of the area 1 m 2 , which is used for referencing and scaling of pictures.

Soil Loss Ratio (SLR)
The sheltering effect of the vegetation is quantified with use of Soil Loss Ratio (SLR). It is determined by referencing the sediment yield from the vegetated plot to the sediment yield from a standard, reference plot maintained as a fallow with seedbed conditions. The reference fallow is cultivated before each simulation according to the procedures stated in the technical standard [26]. First, the topsoil is ploughed into the depth of about 10 cm, aggregates loosened with the rake and levelled with the wooden batten. Large stones and remains of the plants are removed and at the end the surface is carefully rolled with a manual turf roller. Main goal of these procedures is to ensure the surface conditions as close to the reference state as possible for every simulation.

Soil surface roughness and surface changes
The soil surface changes and roughness are evaluated by a photogrammetric technique. The accuracy of derived surface representation depends on the camera specifications. The plausible spatial resolution is about 1-2 mm as reported e.g. in [27] and height accuracy compared to more precise laser scanning may be below 2 mm for smooth surfaces, as reported in [28]. Photogrammetric evaluation is based on the stereo imaging of the field plot surface before and after the rainfall simulation. Since the ground control points (GCPs) are necessary for the image referencing, a dural reference frame with 40 equally distributed GCPs was designed. The frame encloses the maximum area of 1.6 x 1.1 m and is equipped with a mobile cross bar labelled by GCPs which makes the size of the evaluated frame variable. The reference frame is kept in horizontal position with four threaded rods driven into the plot before the rainfall simulation. Sony NEX-5N, interchangeable lens Compact System Camera is utilized for the sensing. The reference frame is removed before the simulation start, whereas the rods are left in the plot to allow the frame being put in the same position after the simulation end. Photographs are always taken with the same camera settings (aperture f-8, exposure time 1/40 and less, focal length 16 or 18 mm depending on selected lens). Thanks to the tarpaulin covering the rainfall simulator the lighting has a diffuse character, which makes the sensing less sensitive to the changes in cloudiness and direct sunlight. The images are taken manually with camera oriented vertically in the height of approximately 1.3 m above the ground and the image distance of 0.3 m. Each captured photo covers the whole frame area and all reference points. Usually more than one stereo pair images are taken for providing larger sample for the evaluation.
Obtained stereo pair images are processed in PhotoModeler Scanner software, which derives accurate, high quality 3D models and enables measurements based on photographs. The process is known as photo-based 3D scanning [29] or image-based modelling [30] and is used for producing a dense point cloud (Dense Surface Modelling) from images of textured surfaces. The multi-sheet camera calibration is executed before the automatic Dense Surface Modelling.
Subsequently, dense surface point clouds created with PhotoModeler Scanner are processed in ArcGIS software. At first the digital surface model (DSM) is generated in the form of triangulated irregular network (TIN) which is then converted to a 1 mm raster DSM. Then the surface raster DSMs created before and after the rainfall simulation are subtracted. This subtraction is used in evaluation of soil volume replacement caused by the rainfall-runoff process. Other analyses are performed on the derived DSMs in order to quantify soil surface roughness parameters such as random roughness, tortuosity, mean upslope depression or other indices. The soil volume replacements as well as the roughness and surface indices are valuable for validation of physically based runoff and erosion models.

SELECTED 2014 STUDY RESULTS
A selection of typical outputs from CTU rainfall simulations is presented. Full results and their in-depth analysis are out of scope of this technical note, provided outputs are to illustrate the most important variables and phenomena discussed in previous sections. During an experimental campaign in June 26 th 2014 the rainfall simulation was carried out on two parallel experimental plots in Bykovicky stream catchment. The field is located in Central Bohemia in the mean altitude of 440 m a. s. l. with mainly well drainable soils -67 % sandy loam, 10 % loamy sand and 8 % sand [31]. The plot steepness is 9 %. The presented experiments were (1) fully grown barley (Hordeum vulgare L.), and (2) reference fallow.

Vegetation cover and soil protection
At the time of the experiments the vegetation has already been fully grown and ready to harvest. Its Canopy Cover was 100%. In this case, therefore, vegetation fully protected the soil surface from direct impact of rain drops.
Measured by LAI-2000 the Leaf Area Index reached 3.8 but as it had been found already in 2013 during comparative measurements, the aging vegetation can affect the results of the LAI-2000 device. This fact is confirmed by the device manual [32] and the published studies [33]. Transmittance of senescent leaves may be larger than that of green leaves, which can result in underestimated LAI [22]. The protective effect of vegetation was assessed by Soil Loss Ratio in 60 minutes from the beginning of the simulation, until then the SLR was close to zero since there was only limited flow of suspended solids recorded in first 40 minutes. The terminal SLR in 60 minutes reached 0,013. Concerning the total average concentrations of suspended solids in runoff from vegetation it reached 0.31 g/l while from the cultivated fallow it was approximately 20 times higher (7.10 g/l).

Suspended solids and phosphorus transport
The total amount of transported suspended solids from fallow plot was 1770 g, approximately 50 times higher than the quantity transported from areas with vegetation (30 g). Both surfaces also varied in terms of the concentration of total and dissolved phosphorus. The concentration of total phosphorus in the cultivated fallow outflow ranged from 2.5 to 7.3 mg/l. In the outflow from the vegetation plot the total phosphorus concentration values were significantly lower (0.4 -1.7 mg/l). On the other hand higher values of the concentration of dissolved phosphorus were measured in the outflow from vegetation (0.136-0.187 mg/l). The concentration of dissolved phosphorus in runoff from the cultivated fallow ranged from 0.071 to 0.088 mg/l.

Soil surface roughness and surface changes
The rectangle of 1.3 x 0.9 m inside the reference frame was analysed for surface changes of the fallow plot. Point clouds (each about a million points) and digital terrain models with a pixel size of 1 x 1 mm were generated before and after the rainfall. The rain and runoff caused a slight decrease in surface random roughness from 17.61 mm to 16.60 mm and the overall surface layer decrease of 2.13 mm as a result of the consolidation of a loosened layers and soil loss.

DISCUSSION AND CONCLUSION
Mobile Rain Simulator represents an important alternative to watching real rainfall-runoff events in-situ [7]. The presented device is a comprehensive system to enable detailed monitoring of a wide range of processes associated with rainfall-runoff episodes. The quest for the maximum amount of the measured quantities is motivated by inter alia, time and financial savings.
Recent experiments also help reveal some shortcomings of the device itself and the procedures applied. In Figure 2 there are evident differences and uncertainties in the pressures in the system. During the experiments, it turned out that the chosen type of the pump during the experiments at its maximum power and any increase in the roughness in pipes is affecting the flow. At the same time, the operator is not able to monitor the minor differences in the pressure in the system, it is therefore necessary to use automatic pressure control device to maintain constant pressure. Because a considerable amount of values is processed, it is convenient to store all the recorded data to by single multichannel data logger for processing. The essential and limiting factor for fully functional use of a similar device is the service staff. The minimum number is a trained and experienced team of five researchers. Another limiting factor is the need for adequate water supply. For rain intensity of 60 mm/hour the consumption is around 1. Last but not least the device is a suitable tool for monitoring the process of transport of nutrients (especially phosphorus) with the water flow. Phosphorus is limiting, on the one hand, the fertility of arable land and the quality of the aquatic ecosystem on the other hand. It is changing the forms from particulate to dissolved, and for this reason, it is very difficult to monitor the major forms of phosphorus during real rainfall-runoff and erosion events. A rainfall simulator is an alternative for these types of experiments under vegetation covers [6], [5] or under alternative erosion control practices [39].
There are many types of simulators. Most of them work with rained area up to 1.5 m 2 , because operation of the DS with larger areas is demanding. The devices are described in detail in a number of studies [40], [41], [42]. The presented equipment thanks to its complexity covers the range of series of published studies that are focused on one or two phenomena related to rainfallflow events. The arrangement of two parallel areas of 16 + 1 m 2 allows monitoring of both rill and interill runoff process simultaneously. The obtained data is a valuable input for the calibration and validation of new and existing mathematical models. Effects of ridge-covering mulches on soil water storage and maize production under simulated rainfall in semiarid regions of China.

ABSTRACT
In recent decades, Fiber Reinforced Polymers (FRP) have shown tremendous potential for retrofitting or repairing existing deficient or damaged concrete structural elements due to their superior properties such as high strength, corrosion resistance and ease of application. However, concern arises about the vulnerability of FRP material to combustion under fire condition, since they are usually applied to the exterior surface of structural members. Damage of the FRP strengthening layer due to high temperature is likely to decrease the load carrying capacity of the columns and threaten the safety of the structure. This paper presents numerical investigation of the behaviour of reinforced concrete (RC) columns strengthened with FRP sheets and insulated by a thermal resisting coating under service load and fire conditions. The finite element numerical modelling and nonlinear analysis are made using the nonlinear finite element analysis software ANSYS 12.1 [1]. The numerical model is verified for several FRP confined and insulated RC columns that have been experimentally tested under service load and standard fire tests in the published literature. The obtained numerical results are in good agreement with the experimental ones regarding the temperature distribution and axial deformation response. Consequently, the presented modelling gives an economic tool to investigate the behaviour of loaded FRP strengthened RC columns under high temperatures occurring in case of fire, if the modelling is verified against experimental works. Furthermore, the model can be used to design thermal protection layers for FRP strengthened RC columns to fulfill fire resistance requirements specified in building codes and standards.

INTRODUCTION
The increasing use of FRP in strengthening applications is due to their high strength, durability and excellent corrosion resistance, in addition to small added thickness and ease of application. However, their inferior performance under fire condition presents a threat to the strengthened member since they are usually applied on the outer surface of the structural elements and the strengthening may be totally lost in case of fire [2]. One of the most popular applications of FRP is strengthening of reinforced concrete (RC) columns. This can be done by two different ways, first by applying FRP sheets to the longitudinal direction of the column in order to provide additional flexural strength, or by applying FRP sheets in circumferential direction in order to provide confining reinforcement which increases both axial load capacity and ductility of the column. Steel spirals or hoops have been widely used in confinement of RC columns, yet they are heavy, difficult to install and corrosive. FRPs provide a competitive solution, being non-corrosive, lightweight and easy to install. Youssef et al. [3] developed a semi-empirical model for FRP confinement for circular and rectangular RC column sections based on large scale experimental program to predict the confined strength ( ′ ) and confined axial strain ( ).
Experimental studies were carried out to investigate the performance of FRP strengthened RC columns under service load and fire condition by several researchers [4][5][6]. In order to provide protection of FRP from fire exposure, a coating layer material of good thermal resisting properties, typically gypsum products, may be placed around the columns. In fire test programs conducted by these researchers, RC columns strengthened by multiple layers of carbon fiber reinforced polymers (CFRP) and protected with coatings of different types and thickness were exposed to standard fire load of ASTM E119 [7] while being subjected to service loads. Experimental and numerical results for both thermal and structural aspects were presented. Using insulation layer with proper thickness was capable of increasing the fire resistance time to over five hours of fire exposure [4].
However, few studies in the published literature addressed numerical modelling to predict the structural or thermal behaviour of FRP-confined RC members subjected to fire with multiple types of protection systems [4,8]. Hence, more investigation work is required to model efficiently the behaviour of FRP-strengthened RC columns under elevated temperatures in order to enable designers to accurately predict the fire resistance time and residual strength after fire exposure and provide a reliable tool for design of thermal insulation layers for these columns.

OBJECTIVE
The present paper aims to study numerically the behaviour of RC columns confined by FRP and thermally protected with insulation material under service load and standard fire test loading. To achieve this goal, numerical modelling by finite elements is made that represents the column geometry, boundary conditions, load variation and considers the variation in thermal and mechanical properties of the different constituent materials with elevated temperature. Numerical modelling and nonlinear analysis are made using the software ANSYS v.12.1.0 [1]. A numerical study is conducted on FRP confined and insulated columns that have been previously tested experimentally under standard fire test. The numerical thermal and structural results are presented and compared to published experimental and numerical results so as to verify the adequacy of the adopted numerical procedure. Finally, the conclusions of the study are given.

VARIATION OF MATERIALS PROPERTIES WITH ELEVATED TEMPERATURE Density
The density ( ) of concrete varies with raising the temperature as given by Eurocode 2 [9], and is shown in Figure 1. The steel reinforcement density is stated by Eurocode 3 [10] to remain constant under elevated temperature. The FRP density stays constant up to approximately 550°C then it has slight decrease. After that, it stays constant until 1000°C (Griffis et al. 1984) [11] as shown in Figure 2. The relations shown in the figures are plotted as normalized density related to density at room temperature.

Thermal Conductivity
The variation of thermal conductivity of concrete with elevated temperature is given by Eurocode 2 [9] and shown in Figure 3. Figure 4 illustrates the variation in thermal conductivity for steel reinforcement with elevated temperatures as given by Eurocode3 [10]. Based on the data provided by Griffis et al. [11], FRP materials have low thermal conductivity and thermal conductivity of FRP seems to decrease with elevated temperatures as shown in Figure 5.

Specific Heat
The variation of specific heat of concrete with elevated temperature is given by Eurocode 2 [9] and shown in Figure 6. Based on the experimental data provided in Eurocode 3 [10], the specific heat of steel with elevated temperature is shown in Figure 7. The complex chemical reactions under high temperature with the composite are the reason behind the irregularity behaviour of the specific heat of FRP material (Griffis et al.,) [11] as shown in Figure 8.  [11] Strength and Modulus of Elasticity The elastic modulus mainly depends on the compressive strength of concrete as given by Eurocode 2 [9] in Equation (1).
As shown in Figure 9, concrete with carbonate aggregate experiences slower degradation rate of its compressive strength at elevated temperatures compared to siliceous aggregate concrete. Variation of tensile strength of concrete with temperature is given by Eurocode 2 [9] as shown in Figure 10. As shown in Figure 11, steel loses more than 50% of its original strength at about 600°C [10]. Likewise, it is clear from Figure 12 that, the modulus of elasticity of the reinforcement steel decreases with elevated temperature.  Figure 13 shows the reduction in strength of CFRP and GFRP with elevated temperature. Likewise, stiffness of FRP materials also has similar rapid degradation with elevated temperature as shown in Figure 14.

Thermal Expansion
The variation of thermal expansion is shown in Figure 15 as given by Eurocode 2 [9] for concrete containing two aggregate types. Thermal expansion of steel reinforcement increases with the increase in temperature as given Eurocode 3 [10] and shown in Figure 16.

Transient Creep
Transient creep is an additional non-recoverable strain that develops only during first time heating of concrete member under compression load as illustrated in Figure 17. Transient creep is load induced thermal strain that happens due to the shrinking of cement matrix and loss of water during heating. In this case total strain has three components, mechanical strain, expansion strain and transient creep strain [13]. It should be noticed that the thermal expansion strain has different direction from the other strains. where is the factor that depends on aggregate type and is equal to 1.8 for siliceous aggregate and 2.35 for carbonate aggregate, is the applied stress, is the compressive strength and is the thermal expansion of aggregate.

Thermal Properties of Insulation Materials
Three types of thermal insulation are used in this numerical study namely; Vermiculite-Gypsum (VG), MBrace and Sikacrete. VG Insulation product is the most favourable material in fire protection manufactured by Tyfo®. This insulation consists of two materials; gypsum and vermiculite. Insulation materials properties variations under elevated temperature are presented in Figures (18-20) [15,16].

Fig. 18 -Variation of density of VG insulation
with temperature [15,16] Fig. 19 -Variation of thermal conductivity of VG insulation with temperature [15,16] The MBrace is a spray applied cementitious mortar based containing various insulating fillers which prevents degradation of concrete mechanical properties above 300 o C. As reported by Chowdhury et al., [5] density and thermal conductivity of MBrace at room temperature were known from manufacturer only and the variation of MBrace thermal properties with elevated temperature was not available, therefore assumptions of the thermal properties are made in this study. The thermal conductivity is assumed to remain constant during duration of fire exposure. Furthermore, since MBrace insulation is cementitious mortar based material, the variation of the specific heat with elevated temperature was assumed to be similar to that of the specific heat of carbonate aggregate concrete as suggested by Lie et al., [17] and shown in Figure 21. The Sikacrete®-213F is a spray applied mortar based protection system with vermiculite and acts as the insulating filler, which gives strong thermal insulation properties. As reported by Masoud et al., [18] the density of Sikacrete is equal to 700 Kg/m 3 and the variations in thermal conductivity and specific heat are given by Eqs. (3) and (4)

MODELLING AND NONLINEAR SOLUTION METHODOLOGY
The general approach for performing the thermal and structural modelling and nonlinear analysis are summarized in Figure 22, the main steps include: 1. Building 3-D nonlinear numerical model of the RC column. This process includes member geometry, appropriate definitions for the materials of concrete, reinforcement steel, FRP and boundary conditions. Two separates finite element models were created, one for thermal analysis and the other for structural analysis.

NUMERICAL VERIFICATION STUDY
A numerical study was conducted for verification of the proposed model for FRP strengthened RC columns subjected to elevated temperatures due to fire. Four studied experimental cases by other researchers are chosen for the verification. Table (1) gives a summary of these cases. The first two columns were tested by Bisby [4], the third column was tested by Chowdhury et al., [5], while the last column was tested by Cree et al., [6].

Element Types for Representation
The element types used to represent the different materials for all columns are given in Table 2. Each type of analysis required specific types of elements in order to get a realistic solution.

Description and Geometry
All columns were experimentally tested in the National Research Council in Canada (NRC) using special furnace which is capable of applying elevated temperature according to standard fire test ASTM E119 [7] and service load via hydraulic jack. All columns have circular section with concrete cover of 40 mm and wrapped with CFRP sheets along their full length. The sections of the columns are symmetric about x and z axes therefore, only quarter of the section was modelled in ANSYS, the finite element meshing is shown in Figure 23.

Materials Properties at Room Temperature
The average concrete compressive strength was taken as given in Table (1), yield strength for longitudinal and transverse reinforcements was 400 MPa for first two columns while they equal to 456 and 396 MPa, respectively for the last two columns. All RC columns were fabricated by carbonate aggregate except for column (3) which was fabricated with siliceous aggregate. The CFRP tensile strength is equal to 1351 MPa for columns (1,2), 3800 MPa for column (3) and 849 MPa for column (4). The stress-strain curve for confined RC section by FRP Youssef et al., [3] was used in this model for different temperature as shown in Figure (24).The rest properties of constituent materials are taken from references [4, 5, 6 and 20].
For RC columns confined with both FRP jacket and steel spiral, the stress-strain curve is influenced by the ratio of lateral confinement pressure of FRP (f1f) to the ratio of lateral confinement pressure of steel spiral (f1s). Approximately, the lateral confinement pressure of confined concrete by both FRP and steel spiral is equal to the sum of lateral confinement pressure of FRP and lateral confinement pressure of steel spiral as proposed by (Lee et al.2004) [19].
The confinement pressure of FRP depends on the strength of FRP sheet which decreases with elevated temperature as discussed earlier. Bisby [4] observed that 90% of the mass of the S-Epoxy was lost in the temperature range of 390°C to 510°C and auto-ignition of the S-Epoxy occurred at approximately 450°C. Therefore, in order to represent this observation in the numerical model, the lateral pressure of FRP was assumed approximately equal to zero after 400°C, whereby the remaining confining pressure mainly comes from steel spirals confinement for higher temperatures.  Figure 25 shows the temperature distribution through the section of column (1) after one, two, three and five hours of exposure. Furthermore, the obtained thermal results using the present model as well as the experimental results obtained by other researchers for the four columns are shown in Figures. 26-29.

THERMAL NUMERICAL RESULTS
For column (1) in Figure 26, the temperature between FRP and concrete changes from room temperature up to 425°C experimentally and 440°C numerically with 3.5% variance after five hours, however, the overall average variance is approximately 24%. As for RFT, temperature increases up to 200°C experimentally and 237°C numerically with 19% variance after five hours, however, the overall average variance is 14%. As for concrete, temperature increases up to 165°C experimentally and 124°C numerically with 25% variance after five hours, however, the overall average variance is 17%. For column (2) in Figure 27, the temperature between FRP and concrete changes from room temperature up to 180°C experimentally and 258°C numerically with 43% variance after five hours, however, the overall average variance exceeds 50%. This happened because the temperature of FRP from first hour to fourth hour remained unchanged below 100°C regardless the continuous heat transfer. As for RFT, temperature increases up to 104°C experimentally and 130°C numerically with 25% variance after five hours, however, the overall average variance is 20%. As for concrete, temperature increases up to 92°C experimentally and 95°C numerically with 3.3% variance after five hours, however, the overall variance is 12%. For column (3) in Figure 28, the temperature between FRP and concrete changes from room temperature up to 480°C experimentally and 460°C numerically with 4% variance after five hours, however, the overall average variance is 12%. As for RFT, temperature increases up to 280°C experimentally and 276°C numerically with 2% variance after five hours, however, the overall average variance is approximately 15%. As for concrete, temperature increases up to 200°C experimentally and 185°C numerically with 8% variance after five hours, however, the overall variance is 27%. For column (4) in Figure 29, the temperature between FRP and concrete changes from room temperature up to 420°C experimentally and 440°C numerically with 5% variance after four hours, however, the overall average variance is 9%. As for concrete, temperature increases up to 180°C experimentally and 220°C numerically with 23% variance after four hours, however, the overall variance is 17%. In general, it can be concluded that the numerically predicted results are in good agreement with experimental ones. The main reason for the variances in prediction temperature mentioned above is generally due to information lack of actual initial values or actual variation with temperature for thermal conductivity, specific heat and density of the used materials in the conducted experiments. These data are not always available from researchers or manufacturers, therefore, standard curves for thermal properties which are discussed in literature have been used in this analysis.

STRUCTURAL NUMERICAL RESULTS
The obtained structural results for the numerically predicted axial deformation under applied service load and fire exposure for specific time are plotted in Figures 30-33 for the four studied columns, and results are compared to the axial deformation measured in published experimental works. From these figures, it can be observed that the numerically predicted results are in good agreement with experimental ones. For column (1) shown in Figure 30, the main deformation is small elongation till first two hours, then elongation strain decreases due to increase of transient creep strain with elevated temperature which may acts against elongation strain until crushing of column. Elongation strain has more influence than transient creep strain in column (2) shown in Figure 31 due to lower exposed of temperature comparing with column (1). As for columns (3,4) shown in Figures 32, 33, columns experienced elongation strain in first two hours then transient creep strain starts to influence the overall deformations.

CONCLUSIONS
The paper presents numerical modelling procedure by finite elements that simulate the performance of RC column strengthened with CFRP sheets and thermally insulated when exposed to standard fire test under service loads. Numerical modelling and nonlinear analysis are performed using the general purpose software ANSYS 12.1. The proposed procedure is verified by comparing the numerical results with experimental results available in the published literature. Based on the obtained numerical results, the following conclusions can be drawn.
1. The numerically predicted results of the proposed model are in good agreement with published experimental results for both thermal and structural aspects. The accuracy of thermal prediction for FRP, RFT and concrete are 78%, 84% and 82%, respectively. Also, for structural failure load prediction, the accuracy is 89%. 2. The proposed model can accurately predict thermal and structural response for different configurations regarding constituent material properties or insulation type and thickness. 3. The proposed model provides a reasonable accurate axial deformation response due to consideration of transient creep strain in the analysis. 4. The axial deformation response for column under fire tend to expand in early stage of fire exposure followed by contraction due to degradation in stiffness and presence of transient creep strain which acts against expansion strain of burned RC column. 5. The proposed model provides an economic tool for check and design of fire insulation layers for FRP strengthened RC columns. It should be noticed that, the modelling needs to be verified against experimental works.

INTRODUCTION
Manned lunar exploration has attracted renewed interest, many countries are preparing their own plans for lunar exploration [1]. The lunar soil will be the first lunar material that explorers and exploration equipment directly contact when they land on the Moon and its mechanical characteristics are important to lunar exploration. Because of the limited amount of lunar regolith brought back to Earth by the Apollo missions, lunar soil simulants, as the substitutions of real lunar soils, was developed to investigate the mechanical properties of real lunar soils [2][3][4][5][6], as well as the interaction behaviour between exploration equipment and lunar soils.
Particle flow code (PFC) numerical simulation as a discrete element method (DEM) is being constantly used in the research of the basic physical and mechanical properties of granular materials. The lunar soil consists of discrete solid particles, which can be classified as a granular material or similar to granular soil [2]. Many researchers investigated the mechanical properties of the lunar soil stimulants using PFC numerical simulation. For example, Hasan and Alshibli [1] carried out PFC three-dimensional (PFC3D) modelling of the tri-axial compression test on lunar soil simulant JSC-1A, and discussed the effects of confining pressure, sample density and environment gravity on the shear strength of lunar soil simulant. Jiang et al. [7,8] used PFC to simulate the biaxial compression test under flexible boundary conditions, and studied the influences of the ground environment (excluding van der Waals force) and lunar environment (including van der Waals force) on the formation process of shear zones. The results indicated that the peak strength and the residual strength of the samples with van der Waals force are relatively high in the compression process, and the lunar environment has a significant impact on the failure mode and characteristics (angle and thickness of shear bands) of the sample. Somrit and Nakagawa [9] simulated the formation process of the cement fusion in the high heat impact zone of the lunar surface when the micrometeorite hit the lunar surface by fitting the impact energy distribution equation into the PFC2D software. Tryana and Masami [10] established the numerical model of uniaxial compression test by considering the parallel cohesion between the particles to carry out the PFC modelling of cement in the lunar soil. The results showed that the cement content has a strong effect on the compression properties of lunar soil. Considering the irregular shape of actual lunar soil particles, Lee et al. [11] described DEM simulations with polyhedral particle shapes to reproduce experimental tests on JSC-1A lunar soil simulant. Katagiri et al. [12] developed threedimensional (3D) grain shape characteristics of returned lunar soil and its numerical simulation by using the image-based DEM. The grain shapes were modelled by clumping 10 spheres in the image-based DEM simulations. Matsushima et al. [13,14] proposed a dynamic optimization algorithm method to simulate the complex shape of the actual lunar soil particles, and carried out PFC modelling of study the influence of the micro-parameters such as the indirect contact stiffness, elastic constant and friction coefficient on the forming process of the angle of repose. Furthermore, the interaction behaviour between exploration equipment and lunar soils was also investigated by DEM numerical simulations [16][17][18][19]. Smith and Peng [17] found that surface roughness significantly influence vehicle mobility and efficiency by DEM modelling of wheel-soil interaction over rough terrain. Zou et al. [18] carried out the plate-sinkage test and the track-shoe test to study the stress sinkage and shear properties of the lunar soil. The results indicated the load that the lunar soil can bear increases with the increase of the bearing area, but the bearing capacity weakens with the increase of the porosity. Li et al. [19] simulated the interaction between lunar regolith and wheel by considering the influence of the lunar low gravity environment. As in the DEM numerical simulation work mentioned above, it can be seen that the mechanical properties of lunar soil and the interaction between lunar soil and wheel (probe) are closely related to the lunar environment and micro-parameters of lunar soil.
The gravity on the Moon is about 1/6 of that on Earth, one of the problems demanding a solution is the mechanical properties of the lunar soil under the microgravity environment [20][21][22][23]. In other words, the lunar soil is subject to very low confining stresses under a microgravity environment conditions which the mechanical properties may differ from those observed under regular tests. In our previous work [24], the mechanical properties and deformation behaviour of QH-E lunar soil simulant at the low confining stresses and those of QH-E at the conventional confining stresses were systematically studied by laboratory test. The research results ascertained the mechanical properties of QH-E lunar soil simulant at lower confining stress and provided a reference and standard for the selection of micro-parameters of the DEM numerical simulation.
In this paper, the tri-axial compressive mechanical properties of QH-E lunar soil simulant at the low confining stress is further investigated by DEM simulation, and compares the simulation results with our previous experimental data [24]. The objective of the present work is to determine the mechanical properties and failure modes of QH-E at low confining stresses, and explore the effects of micro-parameters (porosity, friction coefficient) on mechanical properties. Especially, the failure processes and failure modes of QH-E under different low confining stress conditions are discussed based on DEM simulation.

Selection of particle size and generation of numerical model
In our previous experimental tests, a type of lunar soil simulant named as QH-E was developed by Tsinghua University to investigate the properties of the lunar soil [24]. In the present DEM simulations, the same physical parameters of QH-E lunar soil simulant are selected to further investigate the mechanical properties of QH-E by comparing with our previous experimental results. The grain-size analysis of QH-E was conducted according to ASTM D422-63 for particles greater than 0.075 mm in diameter, and the testing result is shown in Figure 1. As we all know in many DEM simulations, particles are commonly approximated by spherical or ellipsoidal shaped particles to keep the required computational resources manageable, and the spherical particle shape greatly simplifies simulations and accommodates the maximum number of particles for any given central processing unit (CPU) execution time budget [25]. In the present simulation, the shapes of particles are also considered to be spherical with the grain size of 1~1.5 mm for convenience of calculation and easily extract the fundamentals of deformation behaviour and failure modes of QH-E at the low confining stresses.

Fig. 1 -Particle size distribution curves of QH-E sample, the upper boundary, lower boundary, and
average of bulk of Apollo samples [24].
Triaxial compression simulations are performed on the cylindrical soil samples using PFC3D numerical simulation following two steps: (1) inputting command stream in PFC3D to generate a 140 mm height and 70 mm diameter flexible cylinder and two rigid loading plates, which has the same size as the experimental sample. (2) Particle aggregation of simulation model in accordance with the aggregation of lunar soil simulant through PFC3D embedded fish language programming, as shown in Figure 2. In the present simulations, the initial porosity of simulation model is 0.38, and the particle size is the range of 1~1.5 mm random distribution, 40980 particles are produced in the PFC3D model. So large number of particles maybe affect the calculating speed and computational time to a certain extent, but the stress-strain and volumetric strain curves outputted from PFC3D are more uniform.

Contact constitutive model
In the geotechnical problem of lunar soil PFC numerical simulation, the contact constitutive model is generally adopted to calculate the interaction force between particles in the PFC software. The normal and tangential contact stiffness of the particles can be determined according to their own stiffness and contact mechanics model. Due to the fact that lunar soil has a certain amount of apparent cohesion [26,27], parallel-bonded model is used to study the tri-axial compression mechanical properties of lunar soil simulant, which has been successfully applied to simulate mechanical properties of lunar soil simulant [28,29]. Figure 3 shows schematic diagram of parallel bond model and its force and torque distribution. The contact part between the two elements is considered to be a disc, which is expressed in a spring with normal and tangential stiffness. This group of springs are uniformly distributed in the contact plane. Due to the existence of the parallel bond stiffness, after this kind of bond is produced, the contact at the relative motion in the adhesive materials induces a force and a torque, the forces and moments acting on the two bonded particles and related to the maximum normal and tangential stress of the material bonded edges. If either of the maximum stress exceeds the corresponding bond strength, the parallel bond is broken. Parallel bond is defined by five parameters: normal stiffness , tangential stiffness , normal strength , tangential strength and adhesive disc radius , respectively. The normal and tangential stresses of parallelbonded model are expressed as follows: The increment of normal and tangential components of the contact force:

Fig. 3 -Schematic diagram of parallel bond model and its force and torque distribution
(1) where A is the contact area; is the normal and tangential increment for the contact displacement within the time step, respectively. The increment of the contact displacement The contact moment is divided into two parts, the torque and the moment: The torque and moment are described as: where J is the moment of inertia of the disc cross section, I is the moment of inertia of the contact surface along the s i   direction axis and via contact point.

Micro-parameters calibration for PFC simulation
Model parameters required in a PFC numerical simulation can be classified into two categories: 1) physical parameters characterizing the geometrical size of the sample and the experimental conditions; 2) microscopic parameters characterizing particles and their contact mechanics model. The values of the physical parameters can be directly determined by the actual situation of the experiment, while microscopic parameters need to be calibrated which would be a very tedious and time consuming process. The simulation of the tri-axial test (in 3D) or the biaxial test (in 2D) is the basic method to determine the microscopic parameters of the PFC model. Through repeatedly adjusting the input of the microscopic parameters of the model, the simulation results are in agreement with the experiment results as much as possible, so as to determine the values of the parameters, and to further carry out other numerical simulation research. The physical and microscopic parameters of the QH-E lunar soil simulant were calibrated as shown in Table 1. The two other parameters of cohesion and internal friction angle need to be determined by the stress strength values of the experimental results. According to the experimental data draw the Mohr's circle and their common tangent, and then obtain cohesion and internal friction angle. The cohesion is 3.1kPa (c =3.1kPa) and internal friction angle is 51.55 0 (  =51.55 0 ). Figures 5(a)-(c) show the deviatoric stress versus axial strain relationship for QH-E samples at the low confining stresses of 6.25, 12.5, and 25 kPa, respectively. The stress-strain curves can be divided into three stages: hardening stage before the peak deviatoric stress, rapid softening stage after peak deviatoric stress, and the residual stress stage after softening. The peak deviatoric stress increases with the increase of the confining stress, and the amount of postpeak softening also increases with the increase of the confining stress (see Figure 5). It is also shown that QH-E lunar soil simulant displayed a significant strain-softening behaviour associated with the development of a shear band as shown in Section 4. The strain softening contributed to an approximatively equal residual strength, and the residual strength envelopes of QH-E are approximately linear. With the increase of confining stress, the axial strain corresponding to peak stress becomes larger. Comparison of the simulation curves and experimental curves, the deviator stress versus axial strain curves are basically consistent in the elastic deformation stage, the magnitudes of stress in the plastic softening stage have slightly smaller than those of experimental results, but have the same deformation trend as shown in Figure 5. Table 2 shows the values of peak deviatoric stress, residual stress and axial strain corresponding to the peak stress at low confining stresses based on both PFC3D simulation and experimental results. As can be seen in Table 2, the simulated peak deviatoric stress, residual stress and axial strain at peak stress are close to those of experimental results. Moreover, from Figure 5 and Table 2, we found the deviatoric stress-axial strain curves and values of peak stress, residual stress and axial strain are closer to the experimental results at the lower confining stress. It indicated that the proposed parallel-bonded model and calibrated microscopic parameters are able to represent the mechanical behaviour of QH-E lunar soil simulant at the low confining stresses. It is also emphasized that the values of peak deviatoric stress, residual stress and axial strain at peak stress are slightly smaller than those of experimental values, and the deviation gradually becomes larger with the increase of confining stress, which could be related to the model without considering the particle shape features. It is shown that QH-E lunar soil simulant is of angular shape with sharp corners and a highly irregular surface because of the Raymond mill method used to grind QH-E [24]. In the tri-axial compression tests, the irregular particle shape can increase the shear strength [11,30]. Therefore, the simulated values of peak deviatoric stress, residual stress and axial strain at peak stress are slightly smaller than those of experimental results.

Volumetric Strain-Axial Strain curves
Figures 6(a)-(c) show the simulated volumetric strain versus axial strain and comparison with experimental curves at low confining stresses of 6.25, 12.5, and 25kPa, respectively. It can be observed that the larger the confining stress, the smaller the ultimate value of the volumetric strain because of the stronger confining restraint. At low confining stresses of 6.25, 12.5, and 25kPa, the samples displayed highly dilative behaviour since the beginning of loading. The lower the confining stress is, the stronger the dilatancy is, and the greater the volumetric strain is. Further, it can be found that the volumetric strains of the samples at low confining stresses do not tend toward stability, and the samples do not exhibit constant critical state volumes (see Figure 6). The volumetric strains will gradually tend toward stability with the increase of the confining stress. By comparison with experimental curves, it can be seen that the PFC3D simulations of volumetric strain versus axial strain curves show an overall good agreement with experimental curves, and the lower the confining stress is, the simulated curves are closer to the experimental curves. It further reveals that the proposed parallel-bonded model and calibrated microscopic parameters are feasible for studying the mechanical properties of QH-E lunar soil simulant at the low confining stresses.  Figure 7 shows the deviatoric stress and axial strain relationship at different porosities and friction coefficients for the same confining stress of 25.0kPa, respectively. As seen in Figure 7, the peak deviatoric stress and residual stress increased in the stress-strain curves with increasing the friction coefficient, whereas decreased with increasing the porosity, and the initial elastic modulus also increased slightly with increasing the friction coefficient or decreasing the porosity. Zou et al. (2008) studied the effect of micro-parameters on the statics characteristic of lunar soil at the conventional confining stress through PFC3D simulation of tri-axial compression test. Their simulation results show that with increasing the friction coefficient, the peak stress increases obviously, and it increases slightly with increasing the particle contact stiffness. While the peak stress decreases with increasing the porosity. The present results indicate that the effects of the porosity and the friction coefficient on the peak stress of QH-E lunar soil simulant at the low confining stress are the same as those at the conventional confining stress by Zou et al. [29]. Table 3 presents a quantitative relationship for the effects of friction coefficient and porosity on the peak stress, residual stress and axial strain corresponding to the peak stress. The peak stress, residual stress and axial strain corresponding to the peak stress are obviously increased with the decrease of the porosity, the peak deviator stress increased obviously from 87.5 to 466.2kPa, the residual stress increased from 36.5 to 156.5kPa, and the axial strain corresponding to this peak stress also increased from 1.35% to 2.25% when the porosity decreased from 0.42 to 0.3 at the same confining stress of 25.0kPa. The effect of friction coefficient is not as obvious as those of the porosity, the peak deviator stress increased from 172.4 to 205.0kPa and the residual stress increased from 67.5 to 77.6kPa when the friction coefficient increased from 0.20 to 0.40, while the change of the axial strain corresponding to the peak stress is very small (from 1.77% to 1.85%) as shown in Table 3.

Effects of porosity and friction coefficient on stress-strain curves
It can be seen from Figure 7 that the porosity and friction coefficient can obviously affect the peak strength of QH-E lunar soil simulant, and porosity has a stronger effect on the peak strength for tri-axial compression simulation at low confining stress. It indicates that the requirements of macro-mechanical strength can be realized by adjusting the micro-parameters. In general, the initial porosity of the sample is given, so the requirement of the peak strength of tested sample can be achieved by changing the friction coefficient.

FAILURE PROCESSES AND FAILURE MODES AT LOW CONFINING STRESSES
Based on the discussion and analysis of the deviatoric stress versus axial strain and volumetric strain versus axial strain curves, it can be found that the proposed parallel-bonded model and calibrated microscopic parameters are feasible for studying the mechanical properties of QH-E lunar soil simulant at the low confining stresses. In order to more clearly observe the failure processes and failure mechanisms of specimen interior at different low confining stresses, a two-dimensional model is chosen to qualitatively simulate the failure mechanisms by using the above parallel-bonded model and microscopic parameters of QH-E lunar soil simulant at the low confining stresses, and the failure modes are recorded during the whole process of simulation. Figure 8(a) shows the failure process and failure mode at the confining stress of 6.25kPa. A local damage occurs in the upper part of the sample at the axial strain of 0.2%, with the increase of the strain, the damage degree increases gradually, and the failure of QH-E lunar soil simulant is from the edge to the internal zone of the sample. Then, the V-type shear zone (double shear bands) appears in the upper part of the sample at the axial strain of 1.45%, and the V-type shear zone becomes larger with the increase of the strain (at the axial strain of 1.70%) as shown in Figure 8(a). The failure process and failure mode at the confining stress of 12.5kPa have similar results, the V-type shear zone is formed from the edge to the internal zone of the sample, and the only difference is that the V-type shear zone moves a short distance to the right edge of the sample, as shown in Figure 8(b). When the confining stress is 25.0kPa, the V-type shear zone is still observed but move further to the right edge of the sample (see Figure 8(c)). With the further increase of the confining pressure, the V-type shear zone gradually changes to a single shear band at the confining stress of 50kPa, and the single shear band with about 52 o along the horizontal direction (see Figure 8(d)). It reveals that V-type shear zone is the main failure mode at the low confining stress. Then the V-type shear zone gradually changes to a single shear band (along the horizontal direction about 52 o ) with increasing the confining stress. The single shear band with about 52 o along the horizontal direction is the main failure mode at the conventional confining stress. Figure 9 compares the simulated failure modes of QH-E lunar soil simulant at different confining stress with the experimental observations. As seen in Figure 9, the failure modes of QH-E by using numerical simulations are in agreement with the experimental observations. These evolution mechanisms of shear band also successfully explain the reason why QH-E lunar soil simulant exhibits a significant strain-softening behaviour in the deviatoric stress versus axial strain curves.
(a) Failure process at the confining stress of 6.25kPa.
(b) Failure process at the confining stress of 12.5kPa.    At low confining stress, the values of peak deviatoric stress, residual stress and axial strain corresponding to the peak stress show an overall good agreement with experimental results, and these values are closer to the experimental values at the lower confining stress. In addition, the deviator stress vs. axial strain and volumetric stress vs. axial strain curves from the present numerical simulations are also similar to those curves from the experimental tests. Therefore, at low confining stress, the proposed model and calibrated microscopic parameters are suitable for simulation the mechanical properties of QH-E lunar soil simulant.
 Both the porosity and friction coefficient can affect the peak strength of QH-E lunar soil simulant, and porosity has a stronger effect on the peak strength for tri-axial compression simulation at low confining stress. The peak stress, residual stress and axial strain corresponding to the peak stress obviously decrease with the increase of the porosity, while these values slightly increase with the increase of the friction coefficient at low confining stress.  At the low confining stress, the V-type shear zone (double shear bands) is the main failure mode. Then the V-type shear zone gradually changes to a single shear band (along the horizontal direction about 52 o ) with the increase of the confining stress. These evolution mechanisms of shear band are in agreement with experimental observation and explain the reason that QH-E lunar soil simulant exhibits a significant strain-softening behaviour.
The above results based on DEM simulation provide the important information about the strength properties and failure modes of QH-E lunar soil simulant at the low confining stress. However, it is emphasized that the current results are obtained under a special contact constitutive model without considering the shape and size distribution of actual particles. Due to the QH-E lunar soil simulant of angular shape with sharp corners and a highly irregular surface, a more detailed PFC simulation based on a more suitable constitutive model is necessary for a more accurate understanding of mechanical properties of QH-E lunar soil simulant.

INTRODUCTION
Numerous studies were executed in last 20 years to analyse the bending resistance of beam-to-column joints. The component method described in EN1993-1-8:2005 [1] was developed for the beam-to-column steel joints between H and I-profiles and loaded in bending and shear. The resistances are determined for separated components. The compression side resistance of the joint is limited by the shear panel resistance. The influence of the size of the web panel to the quality of prediction by simplified models were studied by Standig [2] and Brandonisio [3]. In reality and in general description of component method the web panel and connections are separated. A simplified procedure is offered in EN 1993-1-8:2005, which is condensing the shear panel and connections into one spring. The shear resistance is described in 6.2.6.1.

EXPERIMENTAL STUDY
Several experiments were examined to obtain the resistance of the column web panel and to verify the joint's resistance. The influence of a high axial force is investigated in [5]. An experimental study of the web panels fabricated from high-strength steel is presented in [6] and [7].
An experimental programme done by Sheer et al. [8] is introduced. A tested specimen from the experimental programme is chosen and used in the numerical study hereafter.

Test details
An eaves moment joint loaded by concentrated force as shown in Figure 1a) was tested to study the behaviour of the column web panel in shear. The results are summarised in [8]. The specimen is made of two welded I section with the dimensions shown in Figure 1b). The height of the section is 306 mm, the column web thickness is 2 mm. The width of flanges is 150 mm and the thickness is 3 mm. The yield stress of the column web fyw = 310 MPa and the yield stress of the beam and column flanges fyr = fyc = 232 MPa were obtained from the tensile test. Fig. 1

VALIDATION
For further investigations, RFEA is developed, see [9] and the design resistance is compared with the experiment. The numerical analyses are completed by the RFEM 5.0 (Dlubal RFEM) [10] finite element program system. The purpose of the validation is to verify the behaviour of the column web panel. A similar study referring to a design of haunches is published by authors in [11].
The materials are endowed with non-linear properties. The ideal plastic model with strain hardening is shown in Figure 2. The geometric details of the model in FEA are taken from the experimental measurement and set according to Figure 1. The welds between plates are neglected in order to simplify the FE model. The numerical model is shown in Figure 3 a). The column is fixed at the end while the beam is loaded by a concentrated force at the free edge. The beam is laterally supported to avoid the lateral torsional buckling.
In the numerical model, 4-node quadrilateral shell elements with nodes at its corners are applied with a maximum side length of 10 mm. Six degrees of freedom are in every node: 3 translations (ux, uy, uz) and 3 rotations (φx, φy, φz). Material and geometric nonlinear analysis with imperfections (GMNIA) is applied. Equivalent geometric imperfections are derived from the first buckling mode and the amplitude is set according to Annex C EN1993-1-5:2006 [12]. Large deformation analysis is used and the Newton-Raphson method for solving systems of equations is chosen. The number of loading steps is set to 50, the convergence criteria for tolerance to 1.0% and the maximum number of iterations to 50. The first buckling mode of the model is shown in Figure 3 b). The design resistance obtained in the numerical model is Mu,Rd = 23,75 kNm. The peak load of the specimen obtained from numerical modelling is validated against the experimental results. The resistance measured in the experiment is Mexp = 25,6 kNm. The difference in the load is 7%, which is less than 10% and is considered as a good result. The stresses in the panel at the design resistance are shown in Figure 4.

VERIFICATION
The design procedure for slender plates is described below. The design procedure is verified on the comparison of CBFEM and RFEM models. The buckling analysis is implemented in the software IDEA Statica Connection [13]. The design procedure for class 4 cross-sections according to reduced stress method is described in Annex B of EN1993-1-5:2006. It allows predicting the post buckling resistance of the joints. Critical buckling modes are determined by material linear and geometric nonlinear analysis. In the first step the minimum load amplifier for the design loads to reach the characteristic value of the resistance of the most critical point coefficient αult,k is obtained. Ultimate limit state is reached by 5 % plastic strain. The critical buckling factor αcr is determined and stands for the load amplifier to reach the elastic critical load under complex stress field. The load amplifiers are related to the non-dimensional plate slenderness, which is determined as follows: The reduction buckling factor ρ is calculated according to EN1993-1-5:2006 Annex B. Conservatively, the lowest value from longitudinal, transverse and shear stress is taken. The verification of the plate is based on the von-Mises yield criterion and the reduced stress method. The buckling resistance is assessed as: where γM1 is the partial safety factor.  Figure 5a) and the first buckling mode is in Figure 5b).  Figure 6, Table 1 and chapter 9.2 in [14].  The results of the extended end plate connection loaded by the bending moment are summarised in Table 2 and for loading by a proportional bending moment and shear in Table 3. The failure modes are included. The bending moment is created by the shear force on a lever arm of 1 m.

Sensitivity to the column web thickness
A parametric study of the column web thickness is presented. A welded cross-section is used for the column and the web panel thickness is changing from 5 to 20 mm. The end-plate and column geometry are summarised in Table 4. The difference from the previous example is highlighted. The results for the extended end plate connection loaded by proportional bending moment and shear are summarised in Table 5. Higher resistance is obtained for a column with a thin web up to 12 mm in CBFEM and in EN1993-1-8:2005. The resistance according to the second draft is underestimated. The distribution of plastic strain in the column web panel with the thickness of 8 mm in CBFEM model is shown in Figure 7. Failure modes: Column web in shear CWP Sh., Beam flange in compression -BF Com.

CONCLUSION
The procedure to design slender plates in structural steel joints is proposed. It is possible to use material nonlinear analysis without imperfections and linear buckling analysis to design slender plates in FEM models with a complex geometry. It is proved that the results of RFEM are in good accordance with the experimental results, therefore, it can be used to predict the actual behaviour of the column web panel in structural steel joints. The proposed design procedure is verified on the RFEM.
Differences of results between prediction of bending resistance by proposal of second draft of EN1993-1-8:2017 and EN1993- 1-8:2005 are in all cases up to 10 %. Higher differences are observed for the combination of bending moment and shear, especially between the second draft and CBFEM model.
The FEA with MITC4 2D shell elements (4 notes, 6 degree of freedom, Mindlin hypotheses) may be expected in this case for the column web panel loaded in shear the most accurate solution. The higher differences in prediction are in case out of a practical application. From the results is clear, that the both analytical models give different results to both sides and oscillate round the FEA results.

INTRODUCTION
The study deals with a decrease in the initial clamping force of a thermoplastic pipe flange connection ( Figure 1) due to creep [1,2] in thermoplastic stub ends made from high density polyethylene (PE-HD). The creep in the PE-HD can be described by the creep modulus [3]. The method uses available long term creep modulus values from the EN 1778 standard [4]. The numerical analysis uses the finite element method (FEM) [5,6]. The study presents a proposed method that takes into account the creep modulus, and applies the method to four different detailed variants. The proposed method considers the exposure time, the stationary temperature and the static load. Another technique that takes into account a creep modulus that is dependent on one temperature, stress field and time was published in [7].

FLANGE CONNECTION
As can be seen from the depiction of the parts of the flange connection [8] in Figure 2, the connection is an assembly in which stub ends are welded to pipes which are connected through  Table 1. The fluid in the pipe system is water at an operating temperature of 20 °C with an internal pressure of 1 MPa.

Creep modulus
Creep modulus values (Table 2) for some thermoplastic materials at the exposure time of 1, 10 and 25 years are given in the standard [4]. The creep modulus stated in [4] is defined as the secant creep modulus determined for a given exposure time for a specific stress level and temperature. The creep modulus values of the PE-HD at the temperature of 20 °C for different stress levels and exposure times obtained from the graphs in [4] are given in Table 2.

Tab. 2 -Creep modulus values of PE-HD for 20 °C for the times stated in [4]
Stress level [MPa] Creep modulus values at selected times

Clamping force
The flanged clamping force is considered according to the EN 12573-4 standard [9]. The standard defines the flanged clamping force under assembly and operating conditions.
The flanged clamping force FA [N] under assembly conditions is defined as: where pO [MPa] is the operating pressure. The flange clamping force under operating conditions (2) corresponds to the product of the operating pressure with the sum of the two areas. The first area is defined by the mean gasket diameter. The other area is half of the gasket area determined in a similar manner as in formula (1)

SIMULATIONS
Flange clamping force relaxation at the selected time points is calculated with the finite element method in the ANSYS software package [6]. The finite element mesh is shown in Figure 1 and Figure 2. The SOLID185 linear brick element is used. Surface to surface contacts are defined between the flange connection parts by a pair of contact elements: CONTA174 and TARGE170. Pretension is applied by a PRETS179 element in the middle of the length of each bolt. At this point a cut perpendicular to the bolt axis is defined. The surfaces created by the cut can be shifted towards each other. The result of this shift is a shortening of the bolt, and therefore pretension loading.

Material
The thermoplastic parts of the flange connection are made from high density polyethylene (PE-HD) [8]. The initial creep modulus value E0 = 800 MPa (0 years) at a temperature of 20 °C is considered to be a constant value for all stress levels. The creep modulus curves in relation to stress level at a temperature of 20 °C for exposure times of 1, 10 and 25 years are shown in Figure  3. The curves pass through the creep moduli from Table 2  The linear elastic model is used for all the thermoplastic parts of the flange connection. The elastic modulus value is considered to be the creep modulus value from the curve (Figure 3) for the mean von Mises stress of a part of geometry. The von Mises stress [1,2] is often used to predict the yielding of ductile materials. Creep modulus is determined based on knowledge of the stress level. However, that is generally not known in advance, so an iterative procedure is used.
The weighted mean of the von Mises stress value for the selected part is determined from the mean von Mises stress values on elements with a weight corresponding to the volume of the deformed elements. Therefore, the weighted mean von Mises stress value depends on the selection of the part of geometry, during which the aim is to divide the structure into a small number of parts of geometry with an approximately constant stress for each of them.

Model variants
The selected four variants of the finite element model of the flange connection, each with a different selection of monitored parts of geometry, are:  Variant Aone part: the stub end and the pipe are together.
 Variant Btwo parts: the stub end and the pipe are separate.
 Variant Cthree parts: the pipe is separate from the stub end, which is divided into two parts.
 Variant Deach element is a separate part.
The four variants of the finite element model of the flange connection are solved in the same way in several load steps. First, pretension load equal to the assembly force FA is applied to the bolts in several load steps, with the criss-cross tightening sequence being used. Subsequently, pressure is applied to the internal surface of the pipe, and in the same load step the pretension load of the bolts is increased to the operating force FO. The ends of the model are defined by two cuts perpendicular to the pipe axis. In the first cut all displacements are set to zero. The other cut is closed with the pipe cap. The pipe cap is loaded by internal pressure as an axial load. The axial load corresponds to the assumption in the formula (2).
First, the shortening of the bolts corresponding to the operating force FO is found for the initial creep modulus value E0. Then the obtained shortening of the bolts is set, and the flange clamping force is monitored while creep moduli ( Figure 3) depending on exposure times and stress levels are decreased.
The dependence of the creep moduli on the stress field is solved iteratively. The first iteration uses the initial estimated stress values and this is repeated until the input and output stress field are sufficiently similar. Geometric nonlinear static analyses are used.

Results
The von Mises stress field at the initial time (0 years) is the same for all variants, see The results after the first year of exposure show that due to the creep in the thermoplastic material, the stresses are redistributed, extreme stresses are reduced and the variants are different from one another. According to expectations the greatest difference between extreme values is found for variants A ( Figure 5) and D ( Figure 8). The results also show that in variants B and C the stress field is discontinuous due to the discontinuity of the creep modulus between the individual parts.
The stress fields at both ends of the model are affected by boundary conditions.
The von Mises stress values and the clamping force values for the exposure times are listed in Table 3. The mean von Mises stresses from all of the solved variants and times are compared in three monitored parts corresponding to variant C (Figure 7).
According to expectations the same results were obtained for all variants for the calculation at 0 years (Table 3) with the initial creep modulus value. The highest mean von Mises stress found for the monitored parts is in the stub ends without the part inside the lap joint flanges. The smallest mean von Mises stress is in the part of the stub ends inside the lap joint flanges. The flange clamping force decreased by approximately 40% from the initial value with the most significant decrease occurring during the first year of exposure. Over the same period of time the creep modulus values decreased to approximately 260 to 170 MPa from the initial value of 800 MPa. The subsequent decrease in the flange clamping force for exposure times of 10 and 25 years is not significant.

Tab. 3 -Mean von Mises stress and flange clamping force in the monitored parts
The results of the mean von Mises stress in the monitored parts after the first year of exposure indicate that variant C is different by less than 2.2% and variants A and B are different by less than 9.3% in comparison with the results for the most detailed variant D. At the same time, the results of the flange clamping force value indicate that variant C is different by less than 1.4% and variants A and B are different by less than 6.3% in comparison with the results for the most detailed variant D.

CONCLUSION
Due to creep in the thermoplastic parts of the flange connection, there is a significant decrease in the flange clamping force during the first year of exposure. After this first year, the initial creep modulus decreased by approximately 75% and the flange clamping force by around 40%. According to numerical analysis results it is necessary to check the flange clamping force during the first year of exposure.
The mean von Mises stress and the flange clamping force at all the selected times except the initial time (0 years) for variants A to C were less than 11.5% different in comparison with the most detailed variant D. The results for variant C, which takes into account the three creep modulus values, are comparable with the results for the most detailed and the most demanding variant D. Variant C is sufficiently precise and is usable in practice.
The calculations were performed for the selected times assuming constant temperature and static load. The proposed method uses a creep modulus that takes into account the exposure time, the stationary temperature field and the static load, and can be generally used for other thermoplastic structures.

ABSTRACT
This work focuses on subjective perception of disruptive sounds with different harmonic frequencies and parts. It verifies the core of subjective disturbance on basis of a listening test. The effort is to clarify if there is more problematic the issue of clearly audible beatswhich is more informative contents issueor predominantly low frequencies of sound. Clearly audible impulses in sound records are meant by the definition beats in this articlefor example drums in the recording. The point of interest is mainly low efficiency of sound insulation against low frequency sound as music with significant beats. The research is restricted only on respondents without hearing difficulties. Based on answers from more than 20 respondents there is conclusion from this paper that the research focus is far from clearance and simplicity. Answers are often very inconsistent and complaints about impossible task were noticed multiple times. In general, however, there is very clear that the problem in disturbance is hidden in low frequencies and with audible beats. The low frequencies are the key disturbance agent in this particular matter.

INTRODUCTION
Insufficient sound insulation and its connection with low-frequency sound sources in building acoustics brings a new problem. In reference to this, there is increasing number of user complaints on the clearly audible beats within a noise from neighborhood. Periodically recurring "beats" in music are mostly a low frequency component. The problem is how to effectively decrease transition of these noise components in lightweight multi-layered structures. These lightweight structures are very different from heavyweight in terms of acoustic behavior.
For a multilayer construction we can consider as an illustrative example a wall made of gypsum board on metal frame where the spaces between which are filled with mineral insulationsee Figure 1.

Fig. 1 -Example of lightweight multi-layered structure
The specificity of this type of structure is the mass-pliability-mass resonance and the standing waves in the air gap between the boards. Due to this behavior, this structure generally has higher sound insulation values than a monolayer structure with an adequate basis weight, but there is a problem with certain frequencies, especially in the low frequency range. In most day-to-day activities, this area of low frequencies, where the sound attenuation is small, does not have much relevance to human perception of comfort. However, there is one very important exception to this rule, which has become more and more frequent in recent years due to the development of available electronics. The exception is quality music listening, which contains a large component of low-frequency sound. The increasing number of home cinemas, high-quality speakers with subwoofers and low-frequency speakers, combined with multi-layered designs, have led to a growing number of complaints from residents.
The basic research question of this work is whether the dominant component of sound disturbance is rather low frequency or clearly audible repeating beats. The listening test used to answer this question is narrowed to a healthy audience without any hearing impairment. The finding should help to design structures more efficiently in terms of sound insulating parameters.

Respondents
Since the work is aimed at healthy individuals without auditory pathologies, the respondents were subjected to an indicative audiology examination prior to testing. This examination was conducted under the supervision of the test administrator and is freely available on the website [1] . The original verification of this online audiogram was done by the author of the article which compared the results of this online test and audiogram with the physician -audiologist. Both audiograms differed only in small sessions and therefore it can be concluded that for the purposes of this work the accuracy of freely available audiology investigation of respondents is sufficient. However, it must be emphasized that it is definitely not a full-fledged substitute for an audiogram from a medical professional.

Recordings
The audio recordings used in the listening tests are digitally synthesized from the sound of the undercover music (played on electric keys that was deliberately chosen neutral and not known to the respondents) and drum strokes at 180 BMP (beats per minute). The underlying music remains the same in all the samples, the effect of the percussion instruments, beats. These are equalized to achieve their main frequencies at 50 Hz 500 Hz and 5000 Hz. The other frequencies were subdued from the original, the frequency then slightly increased. These sounds were then paired and submitted to the respondents.
The underlying sound has been muted by 5 dB in comparison to the original, so that the beats cannot stand out and be audible. This background sound was deliberately chosen to be unknown, and thus did not incite other unsupervised unconscious responses of respondents. It was also advisable that this sound itself did not contain clearly audible and recognizable beats. The length of the sound recordings is 7 seconds, long enough to show and record beats, but so short that the whole sound does not shrink (because it is not very interesting, informally or musically).
Sound pressure levels (SPL) for third-octave bands -For beats at 50 Hz, 500 Hz and 5000 Hz -Left always flat, right with filter A:

Listening tests
Respondents who passed an audiology test were presented with listening tests. On the first page there is an introduction to the test and a short questionnaire on the respondent's nationalage, gender, whether he or she has previously participated in a listening test and the last question whether he or she is aware of any hearing impairment. There are 9 screens in which respondents have the task of assigning an appropriate level of volume to the sound judged to match the default sound -see Figure 6. They have the ability to click on the buttons to amplify and attenuate the sound being played and, of course, to play the default sound. After comparing the volume to the best of consciousness, the respondent moves the "OK" button to another form. The entire test is created in C # programming language in Visual Studio 2015 [3] .The whole questionnaire is conceived as anonymous. The sound can be played repeatedly by clicking on the button, but it does not play in the loop -again to better keep the attention of the respondents. The final test form is dedicated to giving feedback and encouraging you to leave a comment. [4] The schema in which the evaluation form was ordered is outlined in Figure 7. In addition to retrospective assessment (e.g., the relationship between 500 Hz and 5000 Hz with 500 Hz), the evaluation of the sound with itself is taken into question. The judged sound has a different volume set from the beginning, so the 50 Hz and 50 Hz comparison results should not be 0 to be correct. The arrow always goes from the default sound (which has a fixed volume to sounds that can be mitigated / amplified by the respondent.) The arrow number then represents the order of the form in the listening test. Tests took place at different places, depending on the respondents' preference. In any case, however, the test was performed within a acoustically insulated room, without other distractions and without significant psychological subtexts.

RESULTS
According the answers from respondent of listening tests there is possible to conclude several things.
In following table there are shown all the raw valid answers. There were also 4 respondents who do not satisfy the preliminary audiology test so their answers are not processed.

Tab. 1 -Sheet with answers from respondents
From the table is clearly visible that the range of answers is quite broad. The conclusion is that it is highly difficult for people to compare different soundsdespite their similarity. Answers depend on each respondent and his preference or maybe better spoken -immunity to particular component of disturbing sound. Therefore, it is not possible to strictly announce that something is true covering all of population.
In the next table there are results already processed by some means. Answers are averaged and adjusted to the same level (nullification of that +5 sound level point in default sound). Also there is added the analyzed sound from sound analyzer, which is very appreciated for objective comparison.

Tab. 2 -Sheet with processed results
The blue columns show the same soundat this point there is agreement within the answers in respondents themselves and analyzer as well. One decibel is close to the Just Noticeable Difference (JND) for sound level [2] . On the other hand there are significant differences in others columns. In general, however there is clearly visible inclination to evaluate sounds with dominant frequencies of 5000 Hz as the most silent ones and 50 Hz as the loudest one -Comparison 50 Hz and 5000 Hz differs about 18,68 points (respectively 14,05 in case of 5000 Hz and 50 Hz) from real A weighted scale based on analyzer. 1 point in answer sheet responds approximately to 1 dB in sound pressure level based on sound analyzer measurement. The difference between 50 Hz and 500 Hz is insignificant. In exact numbers is in average less than 2,5 points.

CONCLUSION
At the first there must be noticed that the population of respondents is very limited and therefore the statistical significance is not so great. There should be more research and study to confirm and specify these results.
The respondents were relatively accurate in evaluation the same soundscomparison 50 Hz and 50 Hz for example. However there was big difference in comparison of other sounds. Also there were multiple complaints of almost impossible task to compare such different soundsalthough only difference was in frequency of the beats.
Answers had anyway shown that in general there is strong inclination to evaluate beats played in higher frequencies as not so loud as beats played on low frequencies. In other words the hypothesis of low frequencies as the dominant disruptive element is confirmed in this paper. If this conclusion is proven right in further studies and research, it should be comprehend into methodic of design of buildings. Based on these results there is crucial way to increase sound reduction especially on low frequencies in order to achieve better insulation capabilities.

INTRODUCTION
SAP (Super Absorbent Polymers) concrete is a kind of porous material with spherical millimetre-size pores inside the concrete [1][2][3]. The average diameter of the saturated SAP is 5mm. In the process of curing, the spherical saturated SAP beads can slowly dehydrate and separate from the cement paste after shrinkage, and eventually leave closed millimetre-size pores evenly distributed in the concrete. SAP concrete have been used in the island far from the mainland, which can substitute the normal coarse aggregate with saturated SAP and reduce the cost by transportation of the normal aggregate. The stress wave propagation characteristic is an important part of the study of the dynamic mechanical properties, and it is also a basic issue in protective engineering application needs to be proved. When the propagation of stress waves in the medium cause plastic deformation of the material, the original elastic stress wave becomes the plastic stress wave. For plastic stress wave, the wave velocity is determined by two parts: one is the density of dielectric material; the second part is the tangent modulus dσ/dε of the plastic part on the stress-strain curve of dielectric materials under dynamic loading. For hardening decreasing material, the tangent modulus of stress-strain curve decreases with the increase of strain. Thus, the wave velocity of plastic stress decreases with increase of stress. The wave profile is becoming more and more flat and finally the wave dispersion appears which is also called the constitutive dispersion [4].
In the case of the hardening increasing material, the elastic-plastic stress wave will become steeper and steeper in the propagation process, and then the shock wave will be formed. For linearly hardened materials, the wave front remains unchanged during propagation. Usually, tangent modulus of the plastic parts of stress-strain curve for a certain material is smaller than elastic modulus, which will cause the propagation speed of plastic stress wave is slower than elastic stress wave, thus the phenomenon of chasing after unloading will appear. When there is interface, the steep wave front will change to form multiple steps due to the reflection of plastic stress wave on the interface and unloading effect of the elastic unloading stress wave. The time history curve of stress pulse appears obvious waveform dispersion phenomenon due to the stepped wave front, and this kind of wave dispersion is different to constitutive dispersion, which is called interface dispersion [5][6][7].
The SAP concrete has the spatial reticulated shell structures, in which there is much free shell surface and three-dimensional distribution of internal stress. There is complex reflection effect for the tress wave propagating between the free surface inside the concrete. Furthermore, the constitutive model of the SAP concrete exhibits nonlinear properties such as plastic characteristics in whole or partial part of the material. [8][9]. So, it is difficult to theoretically conduct an analytical analysis, especially when partial or whole plastic formation appears and it will be more complex than the elastic stress wave to analyse the reflection effect on the interface for plastic stress wave. Therefore, it is difficult to resolve the propagation behaviour of stress wave in the SAP concrete, and it is necessary to carry out the research through the numerical method of finite element. In this paper, the expanded Drucker-Prager model [10][11] has been used with further introduction of the strain rate effect of the mortar matrix to study the stress wave propagation characteristics in SAP concrete.

Simplification of the SHPB test equipment
The simplified SHPB (Split Hopkinson Pressure Bar) experiment device has been used to study the stress wave propagation characteristics in SAP concrete. The simplified calculation model of SHPB experiment device mainly includes the loading stress wave, specimen, incident bar and transmitted bar with uniform cross-section, as shown in Figure 1.

Fig. 1 -Schematic diagram of the numerical simulation for prorogation of stress wave
The circular cross-section of specimen and pressure bars is simplified to square crosssection. The material model for incident bar and transmitted bar used in the simulation is a linear elastic model. The friction coefficient between the surface of pressure bar and specimen is 0, which means the friction effect of the interface is ignored in the calculation. The incident bar and transmitted bar are mainly to load and record the stress wave propagation through the specimen. It is necessary to have a long specimen, thus, a complete stress wave with at least a wave length will exist in the specimen. The length of specimen is set to be 150mm with comprehensive regarding of the complete stress wave propagating in the specimen and calculation cost.

Constitutive model of the mortar matrix
The numerical research conducted in this paper is mainly based on the experimental data of SHPB test of the SAP concrete with different porosities, which is 10%, 20%, 30% and 40% respectively, and under quasi static loading or dynamic loading with three strain rate levels of 70/s, 100/s and 140/s. The corresponding mechanical performance data of the SAP concrete is shown in Table 1.The uniaxial compressive stress-strain curve under quasi static loading of the mortar matrix is shown in Figure 2 (a), and this typical stress-strain curve can be simplified as shown in Figure 2 (b). In the generic numerical simulation software ABAQUS, the expanded Drucker-Prager model can describe the characteristics of the stress-strain relationship shown in Figure 2, and the expanded Drucker -Prager model can also reflect the hydrostatic pressure correlation of the materials.  The pressure bar used in the numerical model includes three different kinds of materials: mortar matrix of SAP concrete, steel fiber concrete, and corundum concrete. The pressure bar of mortar matrix of SAP concrete is mainly used to study the waveform changes before and after the SAP concrete specimens. The pressure bar of steel fiber concrete and corundum concrete mainly used in protective engineering is to study the influence of multilayer material with SAP concrete middle layer on stress wave propagation. When the material of pressure and specimen are the same, the wave impedance of pressure bar and specimen is equal and it can reduce the reflection of stress wave on the interface, at the same time decrease the stress wave attenuation and dispersion in the pressure bar. The linear elastic model was selected for pressure bar and the density, elastic modulus and Poisson's ratio were in accordance with the mortar matrix. When pressure bar using steel fiber concrete and corundum concrete, the material model was also the same linear elastic model, material density and elastic modulus obtained from the research before. For the convenience of description, the pressure of mortar matrix, steel fiber reinforced concrete and corundum aggregate concrete were named pressure bar 1, pressure bar 2 and pressure bar 3 respectively. From the pressure bar 1 to 3, the impedance between pressure bar and SAP concrete specimen increases. The expanded Drucker -Prager model was used to describe the SAP concrete specimen, and the model parameters are listed in Table 2. The parameters for expanded Drucker -Prager model displayed in Table 2 do not consider the strain rate effect of SAP concrete material itself, which needs to be added to the expanded Drucker -Prager model. The strain rate effect of SAP concrete material itself can be obtained by numerical studying the lateral inertial confinement effect and friction effect in SHPB test of the SAP concrete with different porosities. The formula suggested by CEB [12] can be used to fit the function of the strain rate effect caused by material itself and strain rate, by which the fitting parameters β and γ can be obtained for each porosity. Then, the function of parameter β and γ and porosity can be fitted as:

Tab. 2 -Material properties involved in the numerical model
According to the fitting Equation 1 and Equation 2, it can be deduced that parameters β=0.095, γ=0.268 when the porosity is equal to zero, namely this parameter value of β and γ represent the strain rate effect at high strain rates of mortar matrix caused by the material itself. Thus, the variation of dynamic compressive strength increment caused by materials itself of the mortar matrix and strain rate can be obtained as the equation recommended by CEB: The dynamic mechanical behaviour of SAP concrete can be simulated by introduction of the Equation 3 into expanded Drucker -Prager model. However, the Equation 3 may contain some errors coming from the experiment, data processing and numerical model. The strain rate effect parameters β and γ=0.268 for mortar matrix need to be adjusted, based on the stress-strain curves in SHPB test, through trial calculation in order to accurately obtain the stress wave propagation law for SAP concrete with different porosities. The strain rate effect parameters β and γ=0.268 for mortar matrix after calibration are shown in Table 3.

Incident stress wave
The incident stress wave is a triangle wave with the pulse time 40μsor 80μs, the peak stress 5MPa, 15MPa or 25MPa.The six kinds of incident stress waves are shown in Figure 3. For convenience, the pulse time of incident stress wave is represented by "T" and the peak stress by "F", namely the "T40F15" indicates the incident stress wave with pulse time 40μsand peak stress 15MPa.

Mesh generation
To improve the mesh quality for dynamic calculation, the cross -section shape of the pressure bar and specimen should be simplified: 1) for the pressure bar, the circular cross -section was simplified to square cross -section with original cross -section area, length and volume; 2) the cylinder specimen was simplified to rectangular shape with the same thickness and volume. Trial calculation indicated that the mesh quality around the pores was the main factor that influences the precision of dynamic calculation. To avoid hourglass during calculation, the mesh quality must be guaranteed strictly. The structured mesh technology has been and C3D8R (eight nodes three- dimensional stress reduced integral unit) have been adopted to mesh the incident bar, transmitted bar and specimens. Based on symmetry, this paper only establishes a one-fourth model, which requires symmetrical boundary conditions. The mesh parameters for specimen with different porosities are shown in Table 4.

The transmission wave of the incident stress wave T40F5
The transmission stress wave, through pressure bar 1 and SAP concrete specimens with porosity 10%, 20%, 30% and 40% respectively, of the incident stress wave T40F5 with pulse time 40μs and peak stress 5Mpa is shown in Figure 4. The transmission wave of the incident stress wave T40F15 The transmission stress wave, through pressure bar 2 and SAP concrete specimens with porosity 10%, 20%, 30% and 40% respectively, of the incident stress wave T40F15 with pulse time 40μs and peak stress 15 MPa is shown in Figure 5. The transmission wave of the incident stress wave T40F25 The transmission stress wave, through pressure bar 1 and SAP concrete specimens with porosity 10%, 20%, 30% and 40% respectively, of the incident stress wave T40F25 with pulse time 40μs and peak stress 25 MPa is shown in Figure 6. The transmission wave of the incident stress wave T80F5 The transmission stress wave, through pressure bar 1 and SAP concrete specimens with porosity 10%, 20%, 30% and 40% respectively, of the incident stress wave T80F5 with pulse time 80μs and peak stress 5 MPa is shown in Figure 7. The transmission wave of the incident stress wave T80F15 The transmission stress wave, through pressure bar 3 and SAP concrete specimens with porosity 10%, 20%, 30% and 40% respectively, of the incident stress wave T80F15 with pulse time 80μs and peak stress 15 MPa is shown in Figure 8. The transmission wave of the incident stress wave T80F25 The transmission stress wave, through pressure bar 1 and SAP concrete specimens with porosity 10%, 20%, 30% and 40% respectively, of the incident stress wave T80F25 with pulse time 80μs and peak stress 25 MPa is shown in Figure 9.

DISCUSSION
The influence of porosity on ratio of transmitted peak stress Figure 10 shows the ratio of peak stress of transmitted stress wave and incident stress wave (ratio of transmitted peak stress) through pressure bar 1, bar 2 and bar 3 respectively with the change of porosities of specimens. In Figure 10(a) is shown the influence of SAP concrete on stress wave propagation without the influence of impedance matching of different medium. It can be seen from the Figure 10 (a) that: the ratios of transmitted peak stress are all less than 1; the ratios of transmitted peak stress decreases with increase of porosity; the ratios of transmitted peak stress of the incident stress with short pulse time is significantly smaller than that of the incident stress with long pulse time. The SAP concrete has significant attenuation for stress wave, and the attenuation will increase with rise of porosity, while decrease with increase of pulse time. Figure 10 (b) and Figure 10 (c) display the effect of SAP concrete and interface of pressure bar and specimen on propagation property of stress waves considering the impedance mismatch. It can be seen from Figure 10   The influence of wave impedance ratio on ratio of transmitted peak stress Figure11 shows that the ratios of transmitted peak stress vary with porosity when incident stress wave T40F15 and T80F15propagates into specimen through pressure bar 1, pressure bar 2 and pressure bar 3 respectively. The interface impedance ratio of stress wave λ is 1, 1.56 and 5.38 on the interface of specimen and pressure bar 1, pressure bar 2 and pressure bar 3 respectively.
In Figure11 (a) and Figure 11 (b), it can be seen that the ratios of transmitted peak stress decreases with increase of the interface impedance ratio of stress wave λ, especially for the incident stress wave T80F15 with long pulse time the trend is more obvious. The ratio of transmitted peak stress is less than 0.61 and 0.88 for incident stress wave T40F15 and T80F15 respectively when the interface impedance ratio of stress wave λ is equal to 1. When the interface impedance ratio of stress wave λ rises to 5.38, the ratio of transmitted peak stress is less than o.42 and 0.46 for incident stress wave T40F15 and T80F15 respectively. The multilayer materials with SAP concrete as middle layer can greatly improve the attenuation effect of stress wave combined with attenuation effect by SAP concrete itself and adjustment of wave impedance ratio of different materials. The dispersion effect of SAP concrete on the stress wave As can be seen from Figure 4 to Figure 9, almost all the waveform appear dispersion phenomenon after passing through SAP concrete specimens. The wave dispersion effect includes three aspects: one is that the sharp corners of the triangle wave become smooth and even flat; the second is the rise time of stress waves will be extended; the third effect is that the stress wave oscillates after passing through SAP concrete specimens. The SAP concrete has a large number of free surface of spherical pore. Plastic stress waves can be reflected constantly between the interface and the elastic unloading wave has unloading effect. Thus, multiple steps will occur for the steep loading wave and time history curve of transmitted stress wave will appear obvious dispersion phenomenon because of step wave surface, namely the SAP concrete has significant interface dispersion effect for stress waves.

CONCLUSION
SAP concrete has attenuation function on stress waves. The SAP concrete with higher porosity has greater attenuation effect of stress wave. The attenuation effect will increase with decrease of pulse time of stress wave when given certain porosity for the concrete .The ratio of transmitted peak stress decreases with increase of the interface impedance ratio. The multilayer materials with SAP concrete as middle layer can greatly improve the attenuation effect of stress wave combined with attenuation effect by SAP concrete itself and adjustment of wave impedance ratio of different materials. SAP concrete can cause significant dispersion of stress waves, and the wave dispersion effect includes three aspects: waveform become flat; prolong the rising time; waveform oscillates after passing through SAP concrete specimens. The first two dispersion phenomenon can be attributed to the constitutive dispersion, and the third dispersion phenomenon