3 The method of computation

A non-linear pulsation code (Bridger 1984) based on the description by Christy (1967) has been used to model radial oscillations. Equations of motion and continuity were solved using the explicit finite-difference scheme, whereas the energy equation was treated implicitly with the Crank-Nicholson scheme. Each pulsation model represents a fundamental mode radial pulsator and is characterized by four parameters, $L, M, \mbox{~\em$T_{\rm eff}$ }$ and chemical abundances.

The initial model of the stellar envelope was calculated in hydrostatic equilibrium. The equations were integrated inward from the surface, using approximately 5000 integration points, after choosing an outer-zone mass and a zone mass ratio ($\approx$1.1). The number of mass zones considered in the different hydrodynamic models varies between 70 and 80.

Convection was ignored and chemical composition was assumed to be constant. The inner boundary of the pulsating envelope was considered to be a rigid sphere radiating with constant luminosity. At the surface, the standing wave or total reflection boundary condition was applied ( $P_{\rm tot}=0$). The inner boundary must be at a radius of less than one tenth of the stellar radius and must have a temperature ${\sim}1.\times10^7~~\mbox{K}$.

The static model was perturbed using a radius dependent velocity distribution,

$$U_{i,1}={\rm A}_1\left(\frac{R_{i,1}}{R_{n,1}}\right)^{5} + ... ...,1}}{R_{n,1}}\right)^{7} (~\mbox{$\mbox{km}~\mbox{s}^{-1}$ }).$$ (1)

Here R_i,1 is the initial radius in the ${i}{\rm th}$ zone, R_n,1 is the initial radius at the surface and the amplitude A₁, varies for different models and have typical values in the range 5-20   $\mbox{km}~\mbox{s}^{-1}$. This was done to reduce computing time, rather than assume U_i,1=0.

An artificial viscosity Q_i was introduced as a dissipative pressure following the Richtmyer-Morton method, with the equation described by Stellingwerf (1975).

The adopted viscosity parameters were $C_{\rm Q} = 4.0$, following previous work for helium-rich envelopes (see Saio 1985), and $\rm\alpha_v = 0.1$. However, a deeper analysis of the propagation of shock waves in such media would be interesting, since a different viscosity can modify the shape, and the amplitude, of the velocity and luminosity curves.

For temperatures lower than 10 kK the radiative opacity was calculated using hydrogen-deficient Alexander's (1994) opacities, with a standard solar mixture of metals (see http://web.physics.twsu.edu/alex/pub/G93.txt). OPAL opacity tables (Iglesias & Rogers 1996) were used for temperatures higher than 10 kK. These were initially generated in the OPAL web site (http://www-phys.llnl.gov/V_Div/OPAL/) to have solar mixture of metals with enhanced nitrogen or carbon. Final tables were smoothed and interpolated using adapted OPAL software.

The hydrodynamic models are characterized by their maximum kinetic energy per pulsation cycle, E. For a model with n mass zones, the kinetic energy for each temporal step I, is given by

$$E_I =\frac{1}{2} \sum\limits_{i=1}^{n} \Delta M_i ~ U_i^2.$$ (2)

E is the maximum value of E_I over a pulsation cycle. As a consequence of the initial perturbation, Edecreases initially in some models. In contrast, in the case where U_i,1=0, E increases initially. To ensure that our hydrodynamic calculations had converged, models were continued until E reached a constant value and therefore the pulsation had reached a limiting amplitude. However, some unstable models which did not reach limiting amplitudes were used to find the pulsation periods and to establish limits to the amplitudes of the pulsation.

In general, the distribution of the pulsational work integral for models located within the boundaries of the instability region, has a maximum value deeper in the envelope than for models located in the Cepheid instability strip. This is due to the effect of the iron-group bump, as is shown in Fig. 2. The driving region, where the pulsation work is positive, coincides with the region of maximum luminosity (around zone 40) and with the iron-group opacity peak.

Figure 2: The distributions of pulsational work integral, divided by its maximum value (thicker line W), luminosity (L) and opacity (K) at luminosity maximum, for model 20.