Agustin Spalvier
email: aspalvier@fing.edu.uyjesuseirasjesus.eiras.fernandez@gmail.com

Jesus Eiras

Gonzalo Cetrangolo
email: gonzaloc@fing.edu.uy

Vincent Garnier
email: vincent.garnier@univ-amu.fr

Cédric Payan
email: cedric.payan@univ-amu.fr

Internal Pressure Monitoring in a Post-tensioned Containment Building Using Operational Vibration

Keywords: Nondestructive, Nuclear, Structural health monitoring, Nonlinear, Frequency, Civil engineering, Concrete B

published or not. The documents may come

Introduction

Nuclear power generation requires much attention to safety measurements to contain the nuclear material during service and to finally dispose of its residues. Reinforced and prestressed concrete structures, complemented with other barrier materials, have been the traditional solution to safely contain nuclear material. These critical structures need to follow strict building, testing, and monitoring codes to ensure their correct behavior without imposing a hazard to society and environment. Typical tests required by government regulations and technical codes are based on internal pressure tests, where structural integrity and fluid-tightness are assessed [START_REF] Afcen | RCC-CW Recueil de Conception et de Construction--Civil Works[END_REF][2][START_REF] Courtois | Tests and monitoring for PWR containment : an introduction to AFCEN RCC-CW code[END_REF]. These quality assurance tests take several days, time in which the structure is not operational; thus, planned and unplanned maintenance/testing stops are avoided as much as possible [START_REF] Lundqvist | Instrumentation and monitoring of concrete structures in nuclear power plants[END_REF].

In addition to the regular need for monitoring nuclear concrete structures, there is an increasing concern due to the natural aging of concrete material and the fact that plants are reaching the end of their service lives [START_REF] Lundqvist | Instrumentation and monitoring of concrete structures in nuclear power plants[END_REF][START_REF] Naus | Concrete component aging and its significance relative to life extension of nuclear power plants[END_REF][START_REF] Stevens | Materials Reliability Program: Electric Power Research Institute (EPRI) Review of the Japanese Nuclear Operators' (JNOs') Aging Management Plan for Prolonged Shutdown Periods (MRP-435)[END_REF][START_REF] Niel | La poursuite de fonctionnement des centrales nucléaires au-delà de 40 ans[END_REF]. These have prompted several research projects focused on studying many aspects of structural condition assessment and monitoring in order to extend power plants' service life [START_REF] Garnier | Non Destructive Evaluation for containment monitoring[END_REF][START_REF] Sbartaï | Non-destructive evaluation of concrete damages of containment walls in nuclear power plants[END_REF][START_REF] Garnier | Non-destructive evaluation of containment walls in nuclear power plants[END_REF][START_REF] Clayton | Nondestructive Evaluation of Thick Concrete Using Advanced Signal Processing Techniques[END_REF][START_REF] Scott | Internal inspection of reinforced concrete for nuclear structures using shear wave tomography[END_REF]. The goals are to assess the current state of materials and structures, and to anticipate future issues by developing new monitoring techniques, as robust, efficient, and economical as possible.

Nondestructive testing (NDT) and structural health monitoring (SHM) are testing techniques where the structure under study remains fully unaffected. The collected data can then be used to assess the structure's intrinsic properties, for example, to back-calculate material elastic properties, detect deterioration onset or mechanical stress variations [START_REF] Okayasu | Structural health monitoring system for remote inspection of material failure[END_REF][START_REF] Bui | Evaluation of concrete distributed cracks by ultrasonic travel time shift under an external mechanical perturbation: study of indirect and semi-direct transmission configurations[END_REF][START_REF] Jin | Single-impact nonlinear resonant acoustic spectroscopy for monitoring the progressive alkali-silica reaction in concrete[END_REF][START_REF] Hafiz | Monitoring of stresses in concrete using ultrasonic coda wave comparison technique[END_REF]. This is a key condition, highly valued to assess nuclear concrete structures.

Most of the existing research on NDT and SHM applied to nuclear power plant concrete structures focus on assessing concrete's quality, damage detection, and durability [START_REF] Kurtis | Can we design concrete to survive nuclear environments?[END_REF][START_REF] Andrade | Some principles of service life calculation of reinforcements and in situ corrosion monitoring by sensors in the radioactive waste containers of El Cabril disposal (Spain)[END_REF][START_REF] Duffó | Durability of a reinforced concrete designed for the construction of an intermediate-level radioactive waste disposal facility[END_REF][START_REF] Iliopoulos | Detection and evaluation of cracks in the concrete buffer of the Belgian Nuclear Waste container using combined NDT techniques[END_REF]. There are also multiple reports focusing on general monitoring for long-term mechanical properties characterization, creep, stress determination, and prestressing tendon losses [START_REF] Lundqvist | Instrumentation and monitoring of concrete structures in nuclear power plants[END_REF][START_REF] Maruyama | Strength and Young's modulus change in concrete due to long-term drying and heatingupto90°C[END_REF][START_REF] Trivedi | Assessment of in-situ concrete creep: cylindrical specimen and prototype nuclear containment structure[END_REF][START_REF] Li | Nuclear power plant prestressed concrete containment vessel structure monitoring during integrated leakage rate testing using fiber Bragg grating sensors[END_REF][START_REF] Perry | High stress monitoring of prestressing tendons in nuclear concrete vessels using fibre-optic sensors[END_REF][START_REF] Henault | Truly distributed optical fiber sensors for structural health monitoring: from the telecommunication optical fiber drawling tower to water leakage detection in dikes and concrete structure strain monitoring[END_REF][START_REF] Hermand | Optical fiber for 3D imaging of deformations in concrete containers stacked[END_REF].

The specific objective of this research is to develop an NDT/SHM technique capable of monitoring internal pressure changes and pressure leakage of nuclear concrete containments, externally (no need of accessing into the structure's interior), and without affecting the regular structure's service. Such a technique would be useful as a complementary measurement during regularly planned tightness tests. Moreover, it can become fundamental in case of an accident or for long-term structural monitoring, for example for nuclear waste containment, where no other equipment might be available and the access to the interior is impeded or highly unsafe. Another important application is when the pre-installed (embedded) sensors, such as strain gauges, thermocouples, etc. have malfunctioned [START_REF] Simon | Structural monitoring of prestressed concrete containments of nuclear power plants for ageing management[END_REF]. There are several studies where nuclear concrete structures are subjected to internal pressure and various monitoring techniques are analyzed [START_REF] Li | Nuclear power plant prestressed concrete containment vessel structure monitoring during integrated leakage rate testing using fiber Bragg grating sensors[END_REF][START_REF] Simon | Structural monitoring of prestressed concrete containments of nuclear power plants for ageing management[END_REF][START_REF] Praveechandra | Structural monitoring and response analysis of prestressed concrete primary containment of a nuclear power plant[END_REF][START_REF] Parmar | Over-pressure test on BARCOM pre-stressed concrete containment[END_REF] but none of them purse directly our goal. In particular, Hu et al. [START_REF] Hu | Ultimate analysis of PWR prestressed concrete containment subjected to internal pressure[END_REF] studied the mechanical and material behavior of a post-tensioned nuclear power concrete container under internal pressure using a finite element method (FEM) model. They considered similar affecting factors as ours but their goal, to estimate the structure's ultimate pressure capacity, was different from ours.

To pursue our goal, we have carried out experiments on a mock-up structure: a post-tensioned concrete double-walledcylindrical nuclear containment. The structure was gradually subjected to internal pressure during the experiment in several steps which took 75 h. Throughout that time, we collected acceleration signals with two sensors attached to the outer side of the inner structure. With these, we measured the dynamic response of the structure under different pressure levels. This type of testing technique falls within the field of "output-only dynamic analysis" or "operational modal analysis", where the structure's dynamic features are computed from the natural vibrations occurring while the structure is operational; i.e. there is no specific dynamic input prompting the vibration phenomena. Output-only dynamic analysis is a growing field of study [START_REF] Brincker | Introduction to Operational Modal Analysis[END_REF][START_REF] Rainieri | Operational Modal Analysis of Civil Engineering Structures[END_REF][START_REF] Schanke | Operational Modal Analysis of Large Bridges[END_REF][START_REF] Choi | Modal parameter identification of a containment using ambient vibration measurements[END_REF] whose main advantage is its simpler experimental implementation, which does not require the use of bulky vibrational apparatus, it is more economical and does not affect the structure's regular operations. Choi et al. [START_REF] Choi | Modal parameter identification of a containment using ambient vibration measurements[END_REF] carried out experiments in an operating nuclear containment structure based on operational modal analysis techniques. They were able to identify frequencies of vibration and extract material/structural properties, but they did not focus on monitoring pressure or mechanical stresses variations. This research's main originality is the monitoring of a nuclear containment decennial test and not the method per se. An easy testing setup is successfully employed, based on an output-only experimental system, which does not require bulky instruments nor to embed sensors within the structure. This is the first study focused on studying the relationship between modal frequencies of vibration and internal pressure of a large-size nuclear-containing structure, using outputonly methods. This research demonstrates the existence of clear positive correlations between the structure's natural frequencies of vibration and internal pressure. These are promising findings towards developing a robust NDE/SHM technique capable of monitoring internal pressure variation in complex nuclear concrete containments.

Materials and Methods

Vercors: The Reinforced Concrete Container

Vercors inner structure consists of a post-tensioned reinforced concrete cylindrical structure of 15 m diameter, capped with a dome-shaped lid. Walls are 0.40 m thick, and 30 m tall. The structure is founded on a concrete floor 0.70 m thick. Vercors mock-up is a physical model at 1/3 scale of a real nuclear confinement building. An outer independent structure protects Vercors mock-up from external effects, such as wind and drastic temperature variations. Figure 1 presents a photograph of Vercors mock-up from the outside.

Concrete quality control was carried out during its construction using standard cured companion cylinders. These were tested on day 28 after casting and provided results of density, compressive strength [START_REF] Michel | Time-frequency analysis of small frequency variations in civil engineering structures under weak and strong motions using a reassignment method[END_REF], and modulus of elasticity [36]. Table 1 contains the material characterization results, corresponding to the average properties of each concrete batch.

Vercors mock-up comprised a total of 250 t of reinforcing steel, 50 t of post-tensioning steel, and 2500 m 3 of concrete.

Experimental Procedure

The experiment consisted of gradually injecting/releasing air into the containment to increase/decrease its internal pres- Average density ρ 2201 kg/m 3 a f ck,28 is the "characteristic" strength which corresponds to the 5th percentile of the strength data sure. Air pressure and temperature were measured inside the container. The protocol involved six stages, starting at stage 1 at ambient pressure (0 bar). The internal pressure was increased during stages 2 and 4, and decreased during stage 6, at a rate of approximately 0.2 bar/h. Stage 3 began when pressure reached 2 bar; at this point, the pump was turned off for 5 h approximately, until resuming pressure-increase of stage 4. Upon reaching approximately 4 bar, the pump was turned off for 12 h, until starting the pressure-decrease of stage 6, in which the valve was opened. Figure 2 contains the pressure data measured during the entire experiment.

Two accelerometers model PCB 352A24, with nominal sensitivity of 100 mV/g, were utilized during the experiment. Defining the most appropriate number of sensors is a trade-off situation. While using more sensors would improve results, particularly to determine the modal shapes, more costly and cumbersome experimental setups would be required. Because the goal of this research was to focus on the frequencies of vibration and not the shapes, it was decided to maximize the experimental robustness and easiness, so only two sensors were employed. These two accelerometers were attached to the outer face of the structure, both at level 7.5 m and at angles 0°and 36°from the structure's angular reference (the bigger circular through-wall hole-hatch area-is centered at ~265°). Acceleration signals were acquired every 5 min. Each of these time-domain

Signal Processing

Signal processing was carried out in four steps, signal validation, pre-processing, processing, and post-processing.

Signal Validation

Signal validation consisted of inspecting every raw signal in time and frequency domains, to keep only the signals with useful information. As explained below, some stages of the experiment did not yield vibrational responses that could be distinguished from noise. To this end, the total energy of each signal was calculated by numerically integrating the squared signals in time. An "energy threshold" of 8 × 10 -6 V 2 s was selected for accelerometer 1 and 11 × 10 -6 V 2 sfor accelerometer 2. Signals with lower energy than the threshold were discarded during this step. Other signals that exceeded the threshold but had sudden spikes in time-domain, and showed no clear vibrational behavior in frequency domain, were also discarded during this step, as recommended by Brincker and Ventura [START_REF] Brincker | Introduction to Operational Modal Analysis[END_REF].

Pre-processing

Pre-processing consisted in filtering time-domain signals using a Butterworth bandpass filter of 4th order with cutoff frequencies of 3 Hz and 2000 Hz.

Processing

The "processing" step consisted of computing the frequency spectra. The following processing was applied to each timedomain signal obtained with each sensor. First, a 2-s-long, Blackman moving window was applied, with a 75% overlap. Thus, each 5-s-long time-domain signal produced seven 2-s-long time-domain windowed signals. The selection of the type and length of the window was done based on preliminary analyses that supported those features. Each of the windowed signals was padded with 300,000 zeros at their tails. The amplitude-frequency spectrum of each of these seven windowed-padded signals was computed using the FFT algorithm, yielding a frequency resolution of 0.05 Hz. Then, the seven frequency spectra were averaged, obtaining one frequency spectrum per time-domain signal, for each sensor. The corresponding pairs of frequency spectra, associated to each accelerometer were averaged, obtaining one single frequency spectrum associated to every time of signal collection.

Finally, pressure measurements and collected amplitudefrequency spectra were synchronized.

Post-processing

The "post-processing" consisted of finding the frequencies of vibration of the modes of interest. This task was carried out by selecting the peaks within the frequency spectra. The reader should note that there exist more sophisticated methods for frequency tracking [37,[START_REF] Ditommaso | Monitoring the response and the back-radiated energy of a building subjected to ambient vibration and impulsive action: the Falkenhof Tower (Potsdam, Germany)[END_REF]. Nevertheless, given the short duration of the signals (5 s) with respect to the rate of pressure change (0.2 bar/h), the pressure can be considered constant during the signal, so it sufficed using the traditional peak picking from the frequency spectra.

Computation of Frequency vs. Pressure Correlation

Prior to computing the correlations of frequency vs. pressure it was necessary to define which frequencies of vibration should be used. Preliminary analyses detected frequency peaks recurrently observed during the entire test, while some others, were clearly observed in some signals but disappeared in many others. Therefore, we selected frequency peaks at around 30 Hz, 68 Hz, 93 Hz, 105 Hz, 142 Hz, and 192 Hz, to compute the regression, which appear recurrently in most of the signals. Figure 3 shows the frequency spectrum computed for one of the signals collected during stage 6 and the spectrogram of stage 6.

In Fig. 3 we see that there are other peaks besides the ones selected for the analysis. These were either not recurring throughout the experiment, meaning that they did not appear in many signals, or it could be the case of two peaks being very close so their trends could not be successfully distinguished independently, so these were not included to compute the frequency vs. pressure correlations.

Frequency vs. pressure correlations were carried out using the modes corresponding to the mentioned frequency peaks. The final step consisted of computing the correlations between measured frequency peaks and measured pressure.

The procedure for correlation calculation consisted of fitting the data using a multi-variable linear regression, where each nominal frequency mode (f 30 , f 68 , f 93 , f 105 , f 142 , f 192 ) corresponded to each independent variable. An inconvenient of using multi-variable correlation is that the size of the dataset tends to decrease as more independent variables (frequency peaks) are considered in the model, so there are fewer signals where all six peaks have been identified. This issue was overcome by implementing a variant, as explained below. The multi-variable mathematical model is (1)

P a 0 + a 1,30 f 30 + a 1,68 f 68 + a 1,93 f 93 + a 1,105 f 105 + a 1,142 f 142 + a 1,192 f 192 .
This mathematical model would require that in order to fit the data, the six peaks must have been identified. Those signals where one or more peaks had not been identified would not be considered in the model. Thus, the data size used for the fit would become drastically reduced-due to the randomness of the experimental procedure, many of the signals lack of clear identification of one frequency peak-. To overcome this inconvenient, a variant was implemented, where every possible multi-variable linear regression between pressure and any number of existing frequency peaks were computed. Thus, if a certain signal contains all frequency peaks except for one frequency peak, the pressure prediction would be obtained using the multi-variable model that does not include that missing frequency-peak value. To limit the extent of possible results, the developed algorithm was set so that only one or two of the six frequency peaks could be missing. The correlations procedures were all carried out using the frequency-peaks and pressure data obtained during stage 6 (pressure descending).

Once the pressure vs. frequency-peaks correlation was computed, i.e. the regression coefficients of Eq. ( 1)w e r e calculated, we can then use that correlation to find the "predicted pressure" based on the frequency data. The goodness of the correlation was evaluated by calculating the adjusted coefficient of determination R 2 , and the standard error of the regression, S, defined as the average of the distances between the regression's predictions and the actual observations [START_REF] Frost | Resonance frequencies of a spherical aluminum shell subject to static internal pressure[END_REF]. The value S is useful because it has the same units as the predicted magnitude, in this case, bar, and approximately 95% of the predictions are within ± 1.96 × S from the observed magnitude.

Experimental Results and Discussion

During signal validation, it was observed that signals corresponding to stages 1, 3, and 5 contained no useful information and their total energy was lower than the energy threshold. These stages corresponded to time periods where the pump was off (and the valve was closed), where any oscillatory behavior could not be distinguished from noise. Using operational vibrations, the vertical resolution of the acquisition device is set to warrant a good SNR ratio in such a configuration. So outside operational conditions, i.e. during the time periods where the pump was off (and the valve was closed), without the possibility to remotely change the vertical resolution, there was no sufficient vibration amplitude to perform our analyses. From the 1093 acquired signals with each sensor, only around 520 of them (of each sensor) were validated. Those signals not validated were excluded from further analyses.

Frequency-Peak Results with Respect to Internal Pressure

Figure 4 contains the obtained frequencies with respect to measured internal pressures during stages 2, 4 and 6.

The frequency peak results presented in Fig. 4 show consistent behavior during stage 6, depicting a positive correlation between frequency and internal pressure. This trend has been previously observed by Piacsek et al. [40]i na l uminum spherical containers. It should be noted that the total frequency change, between pressure 0 to 4 bar, is very small, in the order of 0.5 to 1 Hz, depending on the mode. The spectral frequency resolution of 0.05 Hz (see Sect. 2.3) allows observation of meaningful trends between internal pressure and frequency for such frequency variations. Results of stages 2 and 4 are not as consistent as in stage 6; some of the frequency peaks show odd behavior, without a clear trend, for example for the modes at ~142 Hz and ~192 Hz (Fig. 4e and at Fig. 4f). However, a hysteretic behavior becomes apparent. We do not know the precise causes leading to such behavior. While this hysteretic behavior could be the result of impaired damage to the structure upon the pressurization, it can be also attributed (in some cases) to the difficulty of tracking the frequency modes during the experiment which frequently appear coupled. Further tests are needed to elucidate the nature of these apparently odd observations. Also, there are several pressure levels in which frequency peaks could not be detected, for example, the frequency peak f 68 (Fig. 4b) was not detected at pressures from 0 to 2 bar, i.e. during stage 2. However, f 30 (Fig. 4a) shows consistent results for all stages 2, 4 and 6. For these reasons, correlation curves were derived from data of stage 6, and specific physical behavior hypotheses are discussed based on peak f 30 results.

Throughout the experiment, the internal temperature varied between 11.5 and 12.1 °C, and these data did not show any correlation with the measured frequencies.

Frequency vs. Pressure Correlation Curves

The multi-variable linear regression coefficients associated to Eq. ( 1), with all six frequency peaks being identified, one frequency peak missing, and two frequency peaks missing, are presented in Tables 2, 3 and 4, respectively.

Figure 5a presents the graph of predicted pressure with respect to measured pressure, where the unity line is super- imposed. There it is observed that the scattered data are close to the unity line, which indicates qualitatively that the prediction is good. This is also sustained by the adjusted coefficient of determination, which was equal to 0.90. The residual analysis is graphically presented in Fig. 5b, which shows no clear trend, depicting another good feature of the correlation. Figure 5c is another representation of the resid- In the histogram of Fig. 5c, we observe that almost all pressure errors are within -1 to 1 bar. The histogram is Gaussian-shaped, with approximately null mean error and a standard deviation of 0.37 bar; around 67% of the predicted pressures have an error between -0.37 and 0.37 bar. The standard error of the regression (S) is another parameter helpful to characterize the goodness of the regression; it tells how "wrong" the model is on average, in units of pressure. This regression has an S value of 0.29 bar. Thus, approxi-mately 95% of the measured pressures fall within ± 0.57 bar (1.96 × S) from the regression line [START_REF] Frost | Resonance frequencies of a spherical aluminum shell subject to static internal pressure[END_REF].

Pressure Leakage Estimation

Experimental results presented in the previous section show a clear correlation between internal pressure and frequency, where frequency increases with increasing internal pressure. This trend is very clear during stage 6 (pressure decreasing) for all six frequencies of vibration (modes) analyzed here, but not as clear during stages 2 and 4 (pressure increasing). Theoretically, we could use this correlation to identify pressure drops due to leakages of the container. In Fig. 2,w e see a clear pressure drop due to leakage during stage 5 when the pump is off. There, pressure drops from approximately 4.25 bar (last pressure reading of stage 4) to 4.05 bar (first pressure reading of stage 6). This 0.20 bar drop is within the variability (0.37 bar of standard deviation) of our correlation computed using all six frequency modes, so we may not expect to obtain an accurate pressure drop estimation. We note that correlation was intended to represent the global frequency vs. pressure relationship throughout the full range of internal pressures. On the other hand, we d os e ei nF i g .4a that the mode f 30 is particularly interesting because it showed higher consistency than the other modes when comparing behaviors of the stages 2/4 and 6. For these reasons we study the possibility of using the mode f 30 to identify internal pressure leakage. Thus, here we focus the analysis on the leakage by comparing the f 30 results at the end of stage 4 against those at the beginning of stage 6. (2) Equation ( 2) can be used to calculate the "predicted pressure" P pred from a measured frequency peak f 30 . Because of the intrinsic variability of the method, we ought to select a set of signals and work with the set's mean value. Now we select a set of nine signals collected during the end of stage 4, associated to pressures higher than 4 bar, and a set of five signals collected at the beginning of stage 6 associated to pressures between 3.90 and 4.00 bar.

Figure 6a shows the measured internal pressure during time, where the superimposed red circles indicate the times in which time-domain signals were acquired. The end of stage 4 and the beginning of stage 6 are indicated with the names: zone S4 and zone S6, respectively. The very end of stage 4, indicated as zone S4-tip, corresponds to a small period of time where the pump starts to turn off; it provided poor results, so those signals were discarded. The mean values of measured pressures in zones S4 and S6 are 4.10 and 4.01 bar, respectively, indicated with green horizontal bars in Fig. 6c.

In Fig. 6b we show the groups of frequency peaks that characterize S4 and S6, respectively. We see that the frequency population of S4 is slightly higher than the frequency population of S6, which is consistent with having a higher pressure in zone S4 than in S6. In Fig. 6cw ep r e s e n tt h e results of predicted pressure P pred obtained by applying Eq. ( 2) to the mentioned frequency data sets. Again, we see an analogous difference between populations as in Fig. 6b. The mean values of the predicted pressures in zones S4 and S6 are 3.91 and 2.87 bar, respectively. These are shown with black bars in Fig. 6c. A t-test of hypothesis on the means of the two predicted pressure distributions were carried out to prove that these mean values are different with a 95% confidence. However, the predicted pressure difference, 1.04 bar, is ten times higher than the actual difference of pressure means, 0.1 bar. It is highly likely that the structural behavior during stage 6, where pressure is decreasing, is different than during stage 4, where pressure is increasing. However, this effect is in fact an indication of pressure leakage so it could actually be used for this purpose.

From this analysis we conclude that the regression made to match the frequency data over a large range of pressures, from 0 to 4 bar, does not have enough resolution to accurately quantify the amount of pressure leakage during stage 5 (with the pump off). To assess the problem of pressure leakage, we should monitor the frequency variations during the leakage per se. However, this was not possible during our experiment because no useful data was acquired during stage 5, while the pump was off. To overcome this problem, we should employ sensors with higher sensitivity or rely on a different dynamic input, so we could follow the frequencies of vibration while leakage is occurring.

Physical Behavior Discussion

In this section, we use the f 30 mode to investigate and discuss possible underlying physical behaviors that may explain the experimental results. The reader should note that this section is a shallow discussion of potential causes explaining our experimental findings, and we do not draw definitive conclusions. For instance, we note that material-parameter or modal-shapes cannot be extracted by fitting our data into a numerical model. The goal of this section is therefore to aid other engineers, other researchers, and even to ourselves, from a "lesson learned" standpoint, to face new problems of structural health monitoring of complex civil infrastructure.

Estimation of Basic Material Parameters

To investigate the experimental findings from a physical standpoint we used a finite element method (FEM) model using the software COMSOL Multiphysics. The first step required defining the basic material parameters, density, Young's modulus, and Poisson's ratio to feed the model. Because the collected data was not enough to invert the elastic problem to allow finding all three basic material parameters, we decided to estimate density and Poisson's ratio using engineering criteria and only invert the problem to find Young's modulus. The structure's mean density was estimated from applying a weighted average of the concrete's density and steel's density, where the amounts of these were known, as mentioned above in Sect. 2. The Poisson's ratio was 0.20 [START_REF] Meseguer | Jiménez Montoya hormigón armado[END_REF]. Young's modulus was found by matching the FEM model frequency of vibration of mode f 30 to the experimental frequency with null internal pressure; in other words, the structure's concrete's Young's modulus was that one in which the FEM model produced a frequency of vibration of the f 30 mode at 30.50 Hz.

Table 5 presents the basic material parameters utilized for the FEM model analysis.

With the available data, we cannot be fully certain that we have selected the correct mode of vibration to match the FEM model and experimental frequencies. However, Lott [START_REF] Lott | Three-dimensional approach of the elastic nonlinear behavior of rocks and concretes[END_REF] studied this same structure using a different technique based on wave propagation phenomena, two years before our experiments, and determined this same modal shape vibrating at 27 Hz, very close to the frequency of vibration result we observed; moreover, the obtained structure's Young's modulus E RC is consistent with the amounts of the constituent materials (steel, concrete and post-tensioning steel), their elastic characteristics, and the structure's age.

Linear and Nonlinear Effects

In the FEM model, the lower face of the concrete floor was restrained from displacing and rotating. The internal pressure was applied onto the inner faces of walls, ceiling, and floor. A linear-elastic material was considered, using the parameters presented in Table 5, and geometric nonlinearity was enabled. Figure 7 contains an image of the FEM model geometry and the resulting mode-shape associated to the f 30 at null internal pressure.

The internal pressure was increased between 0 and 4 bar in steps of 1 bar by doing a parametric study in COMSOL.

The FEM model and experimental results of f 30 frequency variation with respect to internal pressure are presented in Fig. 8.

In Fig. 8a we see that even though the experiment yielded only one distinct frequency mode near 30 Hz, the FEM model produced three modes around that frequency. The frequencies of vibration of these three modes, as well as the experimental frequency, increase with increasing internal pressure. This frequency vs pressure relationship was previously seen by Piacsek et al. [40] in a thin-wall aluminum container; they founded their results based on the effect of geometric nonlinearity: the increase of walls' inplane tension affects the vibration phenomena by increasing the frequency of vibration, such as it occurs in strings under tension. In Fig. 8, the FEM model mode with the highest frequency is associated to a torsional vibration mode, and the other two with a flexural-type mode as depicted in Fig. 7b andc. These latter modes do not yield the same frequency because of the small non-symmetric structure's geometric features. It is noted that as pressure increases the three modes tend to couple, producing a more complex problem in which frequency vs. pressure relationship becomes highly nonlinear during a limited range of pressures. This phenomenon occurs because not all modes are equally affected by pressure increments.

For the inversion problem of the Young's modulus explained in Sect. 4, we decided to match the average of the two flexural modes to the experimental frequency; thus, at zero pressure we see in Fig. 8a the two frequencies of each flexural mode, starting at 30.47 Hz and 30.54 Hz, where the experimental frequency, 30.50 Hz, is their median value. Still analyzing Fig. 8a, as pressure increases both the experimental frequency and the FEM model frequency results increase. The frequency increments in the FEM model, where the material has a linear constitutive relationship, could be associated with the geometric nonlinearity activated by the applied pressure. However, we see that the experimental frequency increases twice as much as the FEM model flexural frequencies. This additional frequency increment cannot be explained, not even by accepting that our assumed linearelastic material parameters of Table 5 are not correct. Some physical causes that could explain this discrepancy are: (1) the existence of undetected modes coupled with the experimental f 30 mode, (2) concrete material behaving elastically but with softening nonlinearity as internal pressure increases, (3) pressure-dependent boundary conditions, or (4) structural stiffening effect due to the compressed air.

Even though we do not have enough information to elucidate which ones of these four hypotheses are actually affecting the results and how much is their influence, we do show that a slight material softening nonlinearity (in addition to the geometric nonlinearity) yields results compatible with the experimental f 30 . Note that this effect opposes the typical material stiffening observed in plain concrete (not post-tensioned) under compression-acoustoelastic effect [START_REF] Spalvier | The stressinduced surface wave velocity variations in concrete[END_REF][START_REF] Chaix | Study of acoustoelasticity behavior of concrete material under uniaxial compression[END_REF][START_REF] Payan | Determination of third order elastic constants in a complex solid applying coda wave interferometry[END_REF][START_REF] Spalvier | Torsional vibration technique for the acoustoelastic characterization of concrete[END_REF][START_REF] Nogueira | Acoustoelastic Response of Concrete Under Uniaxial Compression[END_REF]-which could be explained by stating that microcracks get closed as compression increases, so the material behaves stiffer at higher compression levels (as long as compression remains below ~30% of concrete compressive strength [START_REF] Payan | Determination of third order elastic constants in a complex solid applying coda wave interferometry[END_REF]). If this was the dominant mechanism, increasing internal pressure should decompress concrete and open the microcracks, and therefore reduce its stiffness and its frequency of vibration. But results show the opposite trend: frequency increases as internal pressure increases, implying therefore a softening behavior. The existence of this softening behavior could be justified because of the post-tensioning force, which is compressing the concrete material to stress levels in which slight softening nonlinearity is likely to occur [START_REF] Meseguer | Jiménez Montoya hormigón armado[END_REF][START_REF]Eurocode 2: Design of Concrete Structures-Part 1-1: General rules and rules for buildings[END_REF]. Thus, we ran the same FEM model but considering a Murnaghan material with third-order elastic constants l + 3000 GPa, m + 2000 GPa, n + 1000 GPa. The f 30 frequency vs. pressure results are presented in Fig. 8b. In that figure, we see that the FEM model flexural frequency curves and the experimental curve have a very similar behavior. It should be noted that this apparent consistency does not prove that this is the actual physical cause governing the problem, but it is useful to guide further research. In any event, these results confirm that the geometric nonlinearity alone does not render a full explanation for the observed behavior; nonlinear elastic behavior should also be considered.

Conclusions

This paper shows the feasibility to monitor the internal pressure of a post-tensioned concrete containment mockup (1/3 scale) during a decennial safety review simulation using operational ambient vibrations. It was possible to monitor pressure changes by monitoring the structure's vibration frequency shifts using a reduced number of accelerometers (two) attached to the outer side of the inner wall of the structure and without measuring the dynamic input. This technique falls in the field of the "output-only dynamic analysis" or "operational modal analysis" which is gaining importance for characterization and structural health monitoring of large structures. The frequency peaks associated to different modes of vibration tend to increase as pressure increases and vice versa. A new variant of a multi-variable correlation was computed to fit the frequencies vs. pressure data. This correlation, which had an adjusted coefficient of determination R 2 equal to 0.9 and a standard error of 0.29 bar, allows to robustly predict the internal pressure. A pressure leakage was clearly identified during pressure bearings. However, the correlation computed over a large pressure range (4 bar) was not capable of successfully predicting the absolute pressure leakage (0.1 bar). An FEM model of the structure was built and the frequency vs. pressure was studied focusing on the mode around 30 Hz. The FEM model frequency increased with increasing pressure, probably due to the geometric nonlinearity induced by the internal pressure. However, the experimental frequency results of the mode around 30 Hz increased twice as much as the results yielded by the FEM model with a linear-elastic material. An FEM model including nonlinear material with softening yielded results that are in good agreement with the experimental data, but this may not be the only cause. Other possible causes which may explain this discrepancy include mode coupling, variable boundary conditions, and/or effect of compressed air. Further research is needed to elucidate which ones of these are significant. For future work, we will test higher sensitivity sensors in order to detect structural vibration during pressure bearing stages, we will use special sensors to measure the lower frequency modes and an increased number of sensors to improve identification of the modal shapes.

Fig. 1

 1 Fig. 1 Photograph of Vercors building mock-up

Fig. 2

 2 Fig. 2 Internal pressure over time during the experiment

Fig. 3 a

 3 Fig. 3 a Frequency spectrum of signal collected during stage 6, where arrows indicate the frequency peaks used for the analysis, b spectrogram of stage 6 with spectral amplitude shown in code of colors, in arbitrary units, and arrows indicate the frequency peaks used for the analysis (Color figure online)

Fig. 4

 4 Fig. 4 Frequency-peak results vs. internal pressure results during stages 2 and 4 (blue circles) and stage 6 (red crosses), for nominal frequency peaks a f 30 , b f 68 , c f 93 , d f 105 , e f 142 , f f 192 (Color figure online)

Fig. 5 a

 5 Fig. 5 a Predicted pressure (P pred ) vs. experimentally measured pressure (P exp ) superimposed with the red dashed unity line, b residuals analysis, and c error (residuals) histogram with a continuous green line

Fig. 6 a

 6 Fig. 6 a Pressure measurements in time. The continuous line (in blue) corresponds to the pressure data measured with a manometer, red circles are the interpolated pressures vs. time in which signals were collected, zones S4, S4-tip, and S6 indicate the three groups of signals. b Frequency values of mode f 30 of zones S4 and S6, respectively. c Blue and

Fig. 7 Fig. 8

 78 Fig. 7 Image of the FEM model of Vercors a geometry (with dimensions in meters); and f 30 mode of vibration shape at null internal pressure observed from b front view and c top view

Table 1

 1 Characterization of concrete using standard cured concrete cylinders at day 28 after casting

	Average compressive strength	f cm,28	48 MPa
	Characteristic compressive strength a	f ck,28	42 MPa
	Average static Young's modulus	E cm,28	37 GPa

Table 2

 2 Multi-variable linear regression coefficients between frequency and pressure when all six frequency peaks are clearly distinguished

	Missing	a 0	a 1,30 a 1,68 a 1,93 a 1,105 a 1,142 a 1,192
	f Nom
	None	-451.8 1.050 0.515 1.088 1.078 0.437 0.564
	Table 3 Multi-variable linear regression coefficients between fre-
	quency and pressure when one of the six frequency peaks has not been
	identified
	Missing	a 0	a 1,30 a 1,68 a 1,93 a 1,105 a 1,142 a 1,192
	f Nom
	30	-458.2 -	0.614 1.104 1.183 0.604 0.541
	68	-451.6 1.175 -	1.133 1.178 0.540 0.573
	93	-449.7 1.024 0.543 -	1.334 0.753 0.697
	105	-426.4 1.414 0.701 1.268 -	0.610 0.686
	142	-452.0 1.224 0.612 1.226 1.243 -	0.669
	192	-422.6 1.351 0.590 1.148 1.202 0.761 -

uals (errors between the predicted pressure and measured pressure) expressed as a histogram.

Table 4

 4

	Multi-variable linear regression coefficients between	1st missing f Nom	2nd missing f Nom	a 0	a 1,30	a 1,68	a 1,93	a 1,105	a 1,142	a 1,192
	frequency and pressure when two of the six frequency peaks	30	68	-460.2	-	-	1.174	1.298	0.787	0.537
	have not been identified	30	93	-454.1	-	0.672		1.434	0.935	0.649
		30	105	-436.2	-	0.866	1.311	-	0.823	0.726
		30	142	-463.9	-	0.754	1.339	1.459	-	0.703
		30	192	-434.5	-	0.678	1.13	1.424	0.941	-
		68	93	-449.3	1.185	-	-	1.482	0.849	0.71
		68	105	-422.7	1.69	-	1.384	-	-	0.715
		68	142	-450.9	1.394	-	1.326	1.424	0.745	0.706
		68	192	-420.0	1.464	-	1.200	1.308	0.889	-
		93	105	-419.9	1.465	0.818	-	-	1.012	0.916
		93	142	-449.3	1.274	0.729	-	1.784	-	0.899
		93	192	-414.0	1.571	0.648	-	1.544	1.127	-
		105	142	-420.0	1.730	0.868	1.557	-	-	0.855
		105	192	-386.6	1.795	0.777	1.520	-	0.974	-
		142	192	-411.2	1.786	0.789	1.466	1.585	-	-

Table 5

 5

	Density (ρ RC ), Young's modulus (E RC ) and Poisson's ratio (ν RC ) of the reinforced	ρ RC (kg/m 2 )	(GPa) E RC	ν RC
	concrete material used for the FEM model	2287.5	49.80	0.20

Acknowledgements

We would like to thank Dr. Jean-Marie Henault for his helpful ideas and discussion, to Prof. Alfredo Canelas for his remarks regarding the numerical model, and to Dr. Salvador Villalobos for his insights about nuclear containments' inspection techniques. We would also like to thank the French ANR project "Non-destructive evaluation of containment nuclear plant structures" (ENDE), for funding this work under the RSNR project, and to EDF for making available the measurements on the mock-up "VeRCoRs".

Funding This project was funded by the French ANR project "Nondestructive evaluation of containment nuclear plant structures" (ENDE) within the RSNR project.

Compliance with Ethical Standards

Conflict of interest

The authors declare that they have no conflict of interest.


This paper addresses the formulation and discrete approximation of dynamic contact initial-value problems. The continuous problem is presented in the context of non-linear kinematics. Standard semi-discrete time integrators are introduced and are shown to be unsuccessful in modeling the kinematic constraints imposed on the interacting bodies during persistent contact. A procedure that bypasses the aforementioned difficulties is proposed by means of a novel variational formulation.

Numerical simulations are conducted and the results are reported and discussed.

Introduction

Mechanical contact is encountered whenever two or more bodies physically interact along their boundaries. Contact is of particular interest in numerous engineering appli cations ranging from metal forming and machine design to soil-structure and structure structure interaction under dynamic excitation.

The non-linear character of contact boundary conditions allows for very few inter esting problems to be solved analytically. Since the inception of the finite element method in the late 50s, numerical solutions to contact problems have been intensively investigated by various researchers. Conry and Seireg appear to be the first to treat contact as a quadratic programming problem, [START_REF] Conry | A mathematical programming method for design of elastic bodies in contact[END_REF], while Chan and Tuba first intro duced a general slideline methodology, [START_REF] Chan | A finite element method for contact problems of solid bodies -part I. theory and validation[END_REF]. Moreover, the works of Zienkiewicz on a flexibility approach, [START_REF] Francavilla | A note on numerical computation of elastic contact problems[END_REF], and Hughes on Lagrange multiplier methods, significantly contributed to the development of robust finite element approximations, applicable to the simulation of large-scale problems.

This work is concerned with the time dimension of contact problems. By way of background, continuous non-linear elastodynamics is briefly reviewed in Section 2. Section 3 addresses the two-body contact problem in the continuous setting, while Sec tion 4 presents Lagrange multiplier and penalty formulations in space and semi-discrete integrators in time suitable for finite element approximations. Section 5 investigates standard second-order implicit time integrators of the Newmark family and their per formance in modeling jump conditions in the kinematic fields (e.g., velocities and accel erations) of the interacting bodies. These integrators are shown to produce undesirable oscillatory solutions along the contact surface. The original Lagrange multiplier for mulation of the two-body problem is appropriately extended so that it eliminate this inconsistency. One-and two-dimensional numerical simulations are offered in Section 6 and demonstrate the significance of the proposed formulation. In this section the initial-value problem of non-linear elastodynamics is briefly re viewed. To this end, a deformable body B is identified with an open set n in the R 3 linear space equipped with the standard basis (Ei. E 2 , E3)• It follows that a typical material point of Bis algebraically represented by vector coordinates X = (Xi. X 2 , X 3 ) in the undeformed (reference) configuration. The boundary of the body, /JO, possesses a unique outer normal N at each of its points, provided it is sufficiently smooth (at least pointwise differentiable). Furthermore, assume that there exists an invertible mapping X defined as x = .X(X, t) , where x denotes the position of particle X in the deformed (current) configuration at a generic time t. The displacement field, u, associated with the motion, is introduced according to u(X, t) = X -X and the deformation gradient, F , as

The strong form of the initial-value problem of non-linear elastodynamics is de scribed by the following set of equations: iio(X) on n.

V•P + P o h P o ii on nxRt,
(5)

In the above, uo and u0 denote the initial displacement and velocity fields of the body, respectively, while P is the 1st Piola-Kirchhoff stress tensor. Also, po = P o(X) is the mass density in the reference state and fu, ft are the Dirichlet and Neumann portions of the boundary, on which boundary displacements, ii, and surface tr�tions, T, are specified, respectively. A hyperelastic constitutive assumption is made in the present context. Thus aw P(X , t

) = P o BF , ( 6
)
where W is a strain energy functional, per unit mass, expressed in terms of some invariant measure of deformation. Equations (1-6) constitute Problem (P).

An integral counterpart of Problem (P} with reference to both the temporal and the spatial dimension can be stated by means of Hamilton's law of Varying Action, [START_REF] Hamilton | On a general method in dynamics[END_REF][START_REF] Hamilton | Second essay on a general method in dynamics[END_REF] U is the space of kinematicall y admissible displacements for the given problem given by

U = {u E H1(fl xRt)' u = u on fu X Rt u(X , 0) = Uo , u(X , 0) = iio on n } ,
while ii is the associated space of admissible variations defined as

ii = {ii E n1 (n x Rt)
ii ii(X , 0) = 0 , fI(X , 0) = O on fl} .

After integration by parts on the kinetic energy, Hamilton's law of Varying Action gives rise to a statement of virtual work for the initial-value problem according to

[' { { [poii • ii + poDuW(u) -pob • ii] dV (7) t, lo -fr I T . ii df} dt = 0 v ii E ii , ( t; ' t fl c R+ .
Equation ( 7) constitutes the basis for the space-time discretization of the problem and will be subsequently utilized in the development of finite elements for dynamic contact. The dynamic contact problem distinguishes itself from the static in that inequality constraint conditions hold not only for the displacements along the contacting surfaces, but also for their rates. The impenetrability constraint may be expressed as

[(Xp2 + u(Xp2, t)) -(Xp1 + u(Xp1 , t))] • np1(X, t) = 0, (8)
where np1 = n is the outer unit normal from body B 1 in the current configuration. Time differentiation of (8) results in

[ii(Xp2 , t) u(Xp1 , t)J • np1(X, t) + [(Xp2 + u(Xp2 , t)) -(Xp1 + u(Xp1 , t))J • ilp1(X, t) = 0,
which, taking into account that the unit normal, np1 is, by definition, parallel to the distance between the contacting points, can be readily shown to reduce to

[ii(Xp2, t) -ti(Xp1, t)]•np1(�, t) = 0. (9)
Likewise, a second time differentiation on (8) yields

(ii(Xp2 , t) -ii(Xp1 , t)] • np1 (X , t) = 0 .
(10)

Similar results may be obtained for higher time derivatives of the boundary displace ments. As it has been just illustrated in the continuous setting, equations ( 9) and (10) are obtained directly from the impenetrability condition and, therefore, do not repre sent additional constraints. They have been derived here because of their significance in discrete solutions.

In the case of persistent contact, the two body problem is summarized as follows:

Problem {Pc)

Given state at time t; for 0'', a = 1, 2, solve (1-6) at interval (t; , t1J subject to p g = 0 , p :'.S 0 , g � 0 on C.

Integral formulations of the problem can be obtained by exploiting the statement of Hamilton's law (or, equivalently, the statement of virtual work) presented earlier in this section. A formulation based on Lagrange multipliers reads:

Problem (VL) Find (u1 , u2 , p) E U1 x U2 x P such that f1 1 t { f )pgii". u" + pgDu.oW"(u") -pgb" • u"JdV lt. o=l Jn -f 'i'". u"'dr } dt + t1 f (fi g + pog)dfdt= o �to h� � v (u1 , u 2 , f>) E u1 x u 2 x f:>, (t; , t1J c R+ .
The spaces of admissible and virtual displacements are as in the previous section, while the respective pressure fields is defined as The overall characteristics of any particular Newmark scheme governed by the choice of the above two parameters. In particular, unconditional stability is guaranteed for 2(3 2 'Y 2 0.5, as shown in [START_REF] Wood | Practical Time-stepping Schemes[END_REF]. Furthermore, integration is globally first-order accurate (error O(tl.t 2 )) provided that 'Y = 0.5.

W't = ( -l ) "LP hAfo n , c, w,, = L g(Vh)A(Ji Q = 1, 2.
For the remainder of this section, all subscripts associated with spatial discretization of the various fields w, ill be dropped in favor of those pertaining to the temporal descrip tion. The Lagrange multiplier formulation of the previous section leads to equations of motion written at time tn+I as 2 L[M" a �+1 + Q"(un+i) -W�+I + W't (u �+I ) ] (

) 14
Use of equations (11) and ( 12) in (13) results in elimination of the unknown velocity and acceleration fields, thus yielding a set of equations with the displacement vector at time tn+t as the only unknown, according to 2 .?; L � L � M"u� +1 + Q"(u� +1) + W£ (u� +1) -W� +t

(15) -M"MLl� � u� + Ll� n v� + (1 - 2 2 ,B)�n = 0 .
Equations (14-15) are generally non-linear and have to be solved by an iterative algorithm such as Newton's method or any of its variants.

5

Dynamic contact/release conditions

In the preceding section, Newmark integrators were employed in the solution of the two-body contact problem, where contact/release constraints were imposed on the dis placement fields, in accordance with (8). Unfortunately, in contrast with the continuous case, the velocity and acceleration fields recovered from such an integration scheme do not satisfy the respective rate impenetrability constraints realized by ( 9) and (10). This can be easily seen with the help of a simple example; consider contact between two nodal points, each belonging to one of the contacting bodies, as in Figure 1. The points come into initial contact at time tn and the following conditions hold:

X1 + u� = X 2 + u� = 0 ' v� = Vo v� = 0,
a�= a�= O.

Enforcement of the impenetrability condition on the discrete displacement fields, as

provided that f3 > O. Likewise, subtracting the discrete velocities of the two nodes a.t

t n +t from one another a.nd using a.gain the initial conditions a.t t n gives

( 2 1 . 1 v n +t -V n +d ••n = �t� (ll n +t -llft+t) • n -V o • n,
which, with the a.id of {16), results in

(v � +1 -v! +i> • n = {� -l)v o • n .

{17)

Equations {16) a.nd {17) illustrate the problematic behavior of Newmark integrators, when directly used in simulating dynamic contact/release conditions. It ca.n be readily concluded that the discontinuity in the post-contact velocities is independent of the time step, while the discontinuity in the accelerations is inversely proportional to the time step, thus consistency is not attained as �t -+ 0. Finally, notice that equations {16) a.nd {17) were obtained without any reference to dynamic equilibrium (which, in this case would only determine the exact position of the nodes at time t n +t ).

The above analysis indicates that a special treatment is necessary for an accurate simulation of contact/release conditions. Indeed, the discrepancy in the computed rate quantities, although immaterial within linear elasticity, can be potentially devastating in non-linear problems employing rate-dependent constitutive assumptions. In an earlier attempt, Hughes and coworkers have used the Newmark solution as a predictor to be subsequently followed by a corrector step, in which velocities were matched on the con tact surface by means of a. local wave propagation analysis. Moreover, accelerations were weighted by the (lumped) masses of the interacting nodes, so that dynamic equilibrium be observed after the treatment, see (HTS + 76]. Here a slightly different methodology is proposed based on a priori satisfaction of the impenetrability constraint and its two rate forms, corresponding to equations {8) and {9-10), respectively. To this end, the Lagrange multiplier formulation in Section 3 is augmented by the variational terms

and 2 f '°'{ w "'(v"' -v n "' -�t n {l -"')a"' n Jch� a n +t ' (19) -�t n 1a� +tl • n (i� +1 • n) + A�(-l) "' i� +t • n + Aa(-l) "' a� +t • n}dr = 0,
where w�, w� > 0 are weighting functions. Equations ( 18) and (19} reveal that the Newmark integrator is modified along the cont. act boundary. Particularly, the admissible contact velocity and acceleration fields are, in principle, assumed independent of the displacements, while the Lagrange multipliers Av and Aa enforce the conditions where

[(•)]c = ((•) 2 -(•) 1 ] • n .
The multipliers Av and Aa are viewed as generalized momenta, energy conjugate to the contact velocities and accelerations, respectively. Equation (18) requires that w�[u! +t -u� -�t n v� -�t� (� -f3)a� -�t� f3a� +1] • n -Av = 0 and

w�[ u ! + I -u ! -At n v ! -At ! ( � -,B)a ! -At ! .Ba ! +i J • n + Av = 0
along Ch. Premultiplying each of the above with the weighting function of the other, recalling that (X + U n +iJc = 0, and, finally, differencing, gives

Av = (X + U n )c + At n (v n )c + At � (� -,B)(a n )c I + I wr W?
Similarly, equation (19) dictates that and 0

w�[v � + l -v � -At n (l -1)a � -At n /R � +iJ • n + Aa 0.
Repeating the same process as before yields

I + I wr W1" (20) (21)
In case one of the contacting bodies is rigid, then it is, by assumption, associated with infinite weighting functions, so that (20) and ( 21 L a {[u� + 1 -u� -At n v� -At � (� -,B)a� J • n + (-1)" �; } In the absence of contact, as well as during persistent contact, the proposed treatment reduces naturally to the associated Newmark method. Indeed, in the former case, as meas ( Ch) = 0, all Lagrange multipliers disappear at the outset, while in the latter, equations (20) and (21) require the multipliers to be again identically equal to zero.

The weighting functions w� and w�, are chosen so that the integrands in ( 18) and ( 19) be (virtual) work quantities. Here, they are set to be where m" denotes a characteristic mass quantity of body a.

Whenever release is detected (as defined in [HTS + 76]), equations ( 18) and (19) yield, in direct analogy to initiation of contact

+ Av ( 1 + -) WV w2 v 1 1 + Aa ( 1 + 2) Wa Wa 0 ( 22
) = 0. ( 23
)
In this work, both Lagrange multipliers in ( 22) and ( 23) are set to zero at release, as their appearance was originally identified with the (currently inactive) constraint of equal velocities and accelerations in the direction normal to the contact surface.

Numerical simulations

Numerical solutions for two dynamic contact/impact problems are presented in this section. Computations are conducted within the environment of the fully non-linear general purpose Finite Element Analysis Program (FEAP), briefly documented in Chap ter 15 of (ZT89) and Chapter 16 of (ZT91).

Although implicit time integration schemes are used throughout, an attempt is made to keep the time step as close to optimal, with regard to the explicit stability limit, as possible. Diagonal (lumped) mass matrices are used for all types of elements, as their behavior is considered superior to that of corresponding consistent matrices whenever, as in given class of problems, the associated kinematic fields exhibit discontinuities (or near-discontinuities) in time. Lumped masses are also used as weighting functions of the modified Newmark impact treatment. Bar 1 is discretized uniformly by 100 elements and Bar 2 by 70 elements, so that, for time step At= 0.142857, both are integrated in time optimall y. Plots of the displacements, velocities and contact force at the contacting tip are shown in Figure 3. Excellent agreement with the exact solutions is exhibited for all three fields.

Cylinder on rigid surface

An elastic cylinder in plane strain is dropped with initial velocity v = -1.0� to wards a rigid surface, see Figure 4. No gravity effects are included in the analysis. The cylinder is modeled by two-dimensional 9-node bi-quadratic elements and the mesh is shown in Figure 5. The elastic constants for the cylinder are E""1 = 5. x 10 2 , """' = 0.3 and the mass density is p = 1.0. The Newmark parameters are set to (3 = 0.25 and 'Y = 0.5. The deformed meshes corresponding to selected stages of the analysis conducted with time step At = 0.05, are illustrated in Figure 6. The history of total force developed due to contact is plotted in Figure 7. Three choices of time step are considered, namely At = 0.02, 0.05, 0.15. It is evident that At = 0.05 produces slightly smoother results than the smaller time step. This can be explained by observing that the stability limit is, in an average sense, closer to one for the larger time •step than for the small er. The larger time step yields somewhat inaccurate resultant forces, as it fails to capture the abrupt changes in the response.

  tThis work has been supported by the Lawrence Livermore National Laboratory under contract Nr. W-7405-ENG-48. A Finite Element Method for Dynamic Contact Problems R.L. Taylor and P. Papadopoulos Department of Civil Engineering, University of California at Berkeley, Berkeley, CA94720, USA 1 2 •Non-linear elastodynamics

  on ft x Rt,

�

  as follows: Problem (VJ For any given time interval (t; , t1] CR+ , find u EU such that D o. 1t1 f [ - 2 1 po ii •ii -poW(u)] dV dt + 1t1 f pob • ii dV dt t; lo t, lo + 111 r T . ii df dt + P o ii . ii]1 1 = 0 v ii E u . � k,In the above equation, Du denotes the Gateaux derivative in the direction of ii. Also,

3

  Dynamic contact problemIn this section a local formulation of the dynamic contact problem is furnished. Considerations are restricted to the two-body problem, with generalization to multi body contact omitted for simplicity. Observing the notation of the previous section, the boundary of bodies B 1 and B2 can be uniquely decomposed into three distinct regions according to an" = ar� U arr U c , °' = i, 2 , where C is the common contact surface. Contact is said to occur if meas(C) > 0. Gap and pressure functions can be defined along C, by means of projection of boundary points from one body onto the other, see [Pap91], as g p r� ) x (an2 r�) x (an 2 f�) x Rt f�) x Rt ,.....

  P :'.So } and Integral statements such as Problem (VL) generally form the basis of finite element approximations to problems in elastodynamics. Formulations based on various penalty regularizations (in the classical penalty, perturbed Lagrangian or augmented Lagrangian sense) can also be obtained in a straightforward manner.4Spatial and temporal discretizationFor a given time, the displacement field, u, and its rates are approximated according to where A(/) is a matrix of basis functions for the entire domain, normally composed by standard interpolation functions within each element (no summation is implied on I).

  Further, vectors u(/)> v(/) and a(/) consist of all nodal displacements, velocities and accelerations, respectively. In addition, the pressure field is spatially approximated according to P � Ph = A(' I) P(l) • Moreover, subscript h characterizes spatially discretized fields. It follows that Prob lem (VL) may be recast in the form: where and Find (ul , u� , Ph) E U� x Ul x Ph such that Q" = L:: vAfoP(u(n),

  In the above definitions, summation symbols over the domain or parts of its boundary refer to the usual assembly operation.A Newmark scheme is employed for the time integration of the spatially discretized equations of motion resulting from the Lagrange multiplier or penalty formulation of the contact problem. The family of Newmark integrators, as originally suggested in[START_REF] Newmark | A method of computation in structural dynamics[END_REF], imposes dynamic equilibrium at the two end points of the finite time interval tl.tn = tn+I -tn. The general form of the integrator is (12) where (•)n = (•)(tn)• Parameters (3 and 'Y (called the Newmark pammeters) may vary according to 0 ::; (3 ::; 0.5

t= 0 Figure 1 :

 01 Figure 1: Dynamic contact; failure of Newmark algorithms defined in (11), along the normal, n, to the contact surface at time tn+l• yields

  ) are modified accordingly. Velocity and acceleration fields normal to the contact surface are recovered with the aid of the above constraint equations as and v� + I • n = [v� + At n (l -1)a� + At n /R� +il • n -(-1)" � w� a � + !• n = .B

  +iJc -At n /(a n +iJc 1 1

0 AFigure e :

 0e Figure e: Impact of dissimilar bars

Figure 3 :Figure 4 :Figure 5 :Figure 6 :Figure 7 :

 34567 Figure 3: Impact of dissimilar bars; history of tip displacements, velocities and contact force