Yandong Gu

Ji Pei

Shouqi Yuan

Jinfeng Zhang

Wenjie Wang

Multi

Keywords:

Introduction

The low-specific-speed centrifugal pump, whose specific speed is between 30 and 80 (ns=3.65nQ 0.5 /H 0.75 ), is widely used in the agricultural irrigations, power plants and chemical productions etc. According to the traditional design method, the impeller passage is very narrow and long in order to meet the requirements of low flow rate and high head. However, the low efficiency of low-specific-speed centrifugal pumps is urgent to increase.

The present optimization on the efficiency and head for low-specific-speed centrifugal pump mainly focuses on the empirical formulas and experiments, but it is difficult for the designers to select the coefficients of empirical formulas and time-consuming to conduct experiments [START_REF] Yongyan | Model of enlarged flow design for low specific speed centrifugal pump[END_REF][START_REF] Zhang Yongxue | Hydraulic design and performance analysis of low specific speed centrifugal pump[END_REF][START_REF] Yang | Parameter optimization of non-overload centrifugal pump[END_REF] . The optimization techniques on the improvements of performance of pump are well-developed. Kim et al applied factorial design to optimize the impeller of pump and analyzed the effects of impeller's parameters on efficiency with response surface method [START_REF] Kim | Improvement of Hydrodynamic Performance of a Multiphase Pump Using Design of Experiment Techniques[END_REF][START_REF] Kim | Design techniques to improve the performance of a centrifugal pump using CFD[END_REF] . A multi-objective optimization of required net positive suction head (NPSHr) and the efficiency of pump under design flow rate was presented by using metamodels and optimization algorithm [START_REF] Yang | Improved Attempt of Non-overload Centrifugal Pumps Base d o n Fluent[END_REF][START_REF] Park | Multi-objective Optimization by Using Modified PSO Algorithm for Axial Flow Pump Impeller[M]//Simulation and Modeling Methodologies, Technologies and Applications[END_REF][START_REF] Zhang | Multi-objective optimization of double suction centrifugal pump using Kriging metamodels[END_REF] . Kim et al optimized the efficiency and cavitation of a mixed-flow pump by combining design of experiment with regression analysis [START_REF] Kim | A Numerical Study on the Improvement of Suction Performance and Hydraulic Efficiency for a Mixed-Flow Pump Impeller[END_REF] .Wang et al. coupled the multi-objective genetic algorithm with back propagation neural network to optimize a transonic axial compressor [START_REF] Wang | Multi-objective optimization of turbomachinery using improved NSGA-II and approximation model[END_REF] . Cravero et al applied an optimization design approach, which combines the design of experiment, response surface method and artificial neural networks, to a four-stage axial flow turbine [START_REF] Cravero | Design optimization of a multistage axial turbine using a response surface based strategy[END_REF] . Wang et al. built expressions of efficiencies and head based on the hydraulic loss and optimized the expression with adaptive simulated annealing algorithm [START_REF] Wang | Optimization method for hydraulic performance of centrifugal pump at multi-operation points[END_REF] .

In this study, a multi-objective optimization procedure of low-specific-speed centrifugal pump is proposed to improve the efficiency under design flow rate. The optimization combines Optimal Latin Hypercube Sampling (OLHS), Kriging model, multi-island genetic algorithm with numerical simulation. Finally, the inner flow is deeply analyzed to illustrate the improvement of performance of low-specific-speed centrifugal pump.

Model pump and numerical simulation 1.1 model pump

The low-specific-speed centrifugal pump IS50-32-160 is determined as the research subject. The pump consists of the impeller equipped with six straight blades and a spiral volute with rectangle section, all of them are made of PMMA for PIV measurements to shoot the inner unsteady flow, as shown in Fig. 1. The design flowrate Q is 6.3 m 3 /h, and the head is 8m. The rotating speed is 1450r/min and the specific speed ns is 46.5. The main geometrical parameters of the model pump are showed in table 1.

The computational domains of the low-specific-speed centrifugal pump are showed in Fig. 2. The software CFturbo is used to build the impeller domain. In order to avoid the backflow occurring and extending to the whole inlet pipe and the recirculation at the outlet pipe, the length of the suction is extended to five times as long as the inlet diameter of impeller and the pipe whose length is five times the outlet diameter of volute is added to the volute.

Mesh generation

Comparing with the unstructured mesh, structured mesh can accelerate the speed of convergence. The quality and the number of the structured mesh are very significant before conducting the numerical simulation, especially for the optimization procedure of the pump.

The structured grids of computational domains are generated by software ICEM CFD, as showed in Fig. 3. The grids are refined in the near-wall flow regions and the mesh independence check was performed with several numbers of grids [START_REF] Mao | Analysis of inner flow in low-specific centrifugal pump at conditions near pump[J][END_REF] . The number of total grids is approximately 3.1 million. In the steady simulation, the shear stress transport (SST) k-ω turbulence model is chosen and the multireference frame (MRF) technique is applied. The inlet and outlet boundary conditions are respectively set as total pressure and flow rate. The interface between the rotor and stator is defined as the Frozen rotor. The physical time scale is set to 0.1ω, where ω is the angular velocity of the impeller. Root mean square (RMS) residual values of the momentum and mass are set to be less than 10 -5 . The performance results are obtained after about 600 iterations.

In the transient simulation, the boundary conditions is the same with that of the steady simulation. The interface between the rotor and stator is changed to "transient rotor stator". The accuracy of discretization in space is set as second order accuracy and second-order backward Euler scheme is used for the time discretization. The total calculating time for the pump is 12 revolutions and the every step time is 1.149×10-4s, corresponding to the rotating angle of 1°. Therefore, the results of final 3 revolutions are used to analyze the pressure fluctuation distributions in the pump. The steady simulation is used to calculate the performance of pump and the transient simulation is conducted to compare with the inner flow characteristics of the original pump and optimized one. The numerical simulations were carried out on a Dell Workstation with an Intel Xeon CPU E5-1650 (12 CPUs) and a clock speed of 3.2 GHz.

Experimental validation

The performance predicted by numerical simulation is necessarily validated by the experiments before conducting the optimization process. The performance of the low specific speed pump was test in showed. The pressure in the inlet and outlet is measured by the WT200 pressure transducers, whose precisions are 0.1%. The flow rates are measured by the IFC 300 electromagnetic flowmeter, whose precision is 0.3%. The comparison between the simulated head and experimental one is showed in Fig. 5. The trend of the head predicted is similar to the experimental performance. The prediction deviations of head is 3.3% under the design condition, and that under 0.7Qd is 3.9%. Therefore, numerical simulation used in the optimization process is reliable.

Optimization process

The optimization procedure is illustrated in Fig. 6. At first, the design variables that have great influence on the pump's performance are determined. Secondly, Optimal Latin Hypercube Sampling (OLHS) is used to generate different schemes of impeller. The impellers are designed within the design space, and then the mesh generation and steady simulation are executed.

The databases are formed by the design variables and the objectives composed of the hydraulic efficiency and head. The approximation model is constructed and the best combination of parameters is obtained based on multi-island genetic algorithm through solving the approximation model. The advantage of the optimization technique is that it can reduce the time of finding the best combination preciously and reduce the computational resource.

Two objectives

The optimization objectives are not only to improve the hydraulic efficiency under design condition, but also to ensure the head to meet the requirements.

 

1.0 2 2 max ( , , ) Q f b     (1)

Design variables

According to the design experience and parameters of impeller that have great influence on the efficiency and head [START_REF] Shouqi | The Theory and Design of Low Specific-Speed Centrifugal Pumps[END_REF] , the design parameters that blade outlet width b2、blade outlet angle β2 and wrap angle φ are selected as the input variables, as showed in table 2.

Optimal Latin Hypercube Sampling method

It is important to make sure that the database is enough to construct the approximation model, and the design space formed by design variables is feasible. Optimal Latin hypercube sampling method, which is improved from the Latin hypercube sampling, is one of the most effective way to build the design space. And the advantage of the sampling method is that it improves the uniformity of the Latin hypercube sampling and the great ability of space filling can make better fitting functions between the objectives and design variables. In the optimization process, 20 design points was generated by optimal Latin hypercube sampling method and the corresponding numerical results are shown in table 3. It can be found that the efficiency of the ninth scheme is highest, whose head also meet the requirements.

Kriging model and Multi-island genetic algorithm

In this paper, the approximation model, called Kriging model, is used to obtain the best combination of design variables of impeller. The function of the Kriging model is to build the approximate mathematic expressions between the objectives and design variables. By solving the mathematic expressions with optimization algorithm, the best combination of parameters can be obtained within the design space quickly.

Kriging model consists of regression function and relative function. The regression function can be polynomial fitting of zero order, first order and second order. The relative function is Gaussian function, linear function, exponent function and so on. Supposed that design variables x of n dimensions is [ x1, x2 … xn ] T ,the response values y(x) is [ y1, y2 … yn] T ,so the y(x) can be written as followed:

y (x)=f(x) +ε (2) f (x)=αg(x)+h(x)
(3) where f(x) is Kriging model. ε is random error between the value of Kriging model and response and it follows normal distribution. g(x) is the regression function of Kriging model , α is the coefficient of regression function,h(x) is correlation function and Gaussian function is selected.

Each design point in the Multi-Island Genetic Algorithm (MIGA), like other genetic algorithms, is perceived as an individual with a certain value of fitness based on the value of the objective function and constraint penalty. The main feature of Multi-Island Genetic Algorithm is that each population of individuals is divided into several sub-populations called "islands". Some individuals are then selected from each island and migrated to different islands periodically. This operation is called "migration". The advantage of the multi-island genetic algorithm is to avoid the phenomenon of "premature", which exists in the general genetic algorithms. The main parameters of multi-island genetic algorithm is showed in table 4.

Results and discussions 3.1 Accuracy of approximation model

It is very necessary to assess the accuracy of Kriging model before the optimization, thus the R-square error is used to analyze the difference between the efficiency and head calculated by numerical simulation and those predicted by the Kriging model. The value of R-square is over than 0.9, as Fig 7 .  showed, it means that the Kriging model can have a satisfactory accuracy of prediction. Therefore, the Kriging model can be applied obtain fine mathematic relationship between the objectives and design parameters combined with genetic algorithm.

Comparison of optimal parameters

The best combination of impeller parameters is figured out by solving the approximation model with multi-island genetic algorithm. The optimized design parameters for impeller of the centrifugal pump are predicted by the Kriging model showed in table 5. It can be seen that the blade outlet width is optimized to smaller, the wrap angle is a little lager, whereas the blade outlet angle is almost the same with the original one. The hydraulic efficiency of the optimized impeller increases by 3.2% under the design condition, whose head also meet the requirements. Besides, the efficiency and head obtained predicted by Kriging model have a good agreement with those calculated by CFD, with a relative deviation of 0.06% and 3.4%, respectively.

Comparison of steady simulations

To find reasonable factors for the improvement of performance in design flow rate, the inner flow characteristics of the original pump and optimized one have been compared, as showed in Fig. 8. The flow separation phenomena occur on the leading surface of the blades in the original pump, but those in the optimized one is much better, because the original one based on increasing-capacity design method operates under the partload condition. The vortex region in optimized impeller is smaller than that in original one and the flow is improved, which mainly results in the of the efficiency of the optimized one. But some impact loss arise in the volute of the optimized pump, thus the matching between the impeller and volute still need to improve. [START_REF] Cravero | Design optimization of a multistage axial turbine using a response surface based strategy[END_REF] . A lower pressure standard deviation value indicates smaller pressure fluctuation intensity. The pressure on each grid node can be given as pi=p(x, y, z, t), the pressure fluctuation intensity psdv is defined as follows:

1 1 ( , , , ) 360 N i n n p p x y z T N    (4) 2 1 1 ( , , , ) 1 3 6 0 N sdv n n p pxyz T p N            ( 5
)
where T is the period of one impeller revolution, 1080 N  . The comparisons between the pressure fluctuation intensity distribution in the original impeller and that in optimized one are showed in Fig. 9. The pressure fluctuation intensity in these two schemes increases along the flow passage. The minimum pressure fluctuation intensity occurs in the inlet and the maximum exists near to the trailing edge at the pressure side, and that on the pressure side is larger than that on the suction side. It is clear that the area of large pressure fluctuation intensity decreases in the optimized impeller under design condition.

The pressure fluctuation intensity distribution in the volutes are showed in Fig. 10. The pressure fluctuation intensity is relatively large from the tongue to cross section II and becomes weaker from cross section VII to VIII, and the largest one is around the tongue because of the interactions of the impeller and volute. Comparing the original scheme, the pressure fluctuation in the optimized one is obviously weakened, which is beneficial to improve the operating reliability. The monitor A near the tongue is set to record the pressure fluctuation of the last revolution showed in Fig. 10. The Fig. 11 (a) shows that the pressure fluctuation is regular and composed of 6 peaks and 6 troughs in a rotating cycle. The amplitude of pressure fluctuation coefficient in the original scheme is larger than that in the optimized one. The frequency of pressure fluctuation of these two schemes are determined by the blade passing frequency, and the amplitude of the dominant frequency in the original scheme is larger than that in the optimized one, but the amplitude of the second frequency in the original scheme is smaller .

Conclusions

An optimization approach based on numerical simulation, design of experiment, approximation model and multi-island genetic algorithm is proposed to improve the performance of the low-specific-speed centrifugal pump and the inner flow characteristics are analyzed to illustrate the improvement of performance.

(1) Kriging model is effective to obtain the fine relationship between the objectives and design parameters and predict the efficiency and head with a high accuracy.

(2) The hydraulic efficiency of optimized impeller increases by 3.2% comparing to the original one under design flow rate, and the head also meets the requirements.

(3) The velocity flow fields through this optimization are improved. The unsteady pressure fluctuation intensity in the optimized scheme decreases under design flow rate.

Fig 2 .

 2 Fig 2. Computational domain of centrifugal pump

Fig 3 .

 3 Fig 3. Structured mesh of the centrifugal pump 1.3 Numerical simulation Taking a three-dimensional and viscous flow in the centrifugal pump into consideration, 3D Reynolds-averaged Navier-Stokes equations are solved by commercial code ANSYS CFX 14.5.In the steady simulation, the shear stress transport (SST) k-ω turbulence model is chosen and the multireference frame (MRF) technique is applied. The inlet and outlet boundary conditions are respectively set as total pressure and flow rate. The interface between the rotor and stator is defined as the Frozen rotor. The physical time scale is set to 0.1ω, where ω is the angular velocity of the impeller. Root mean square (RMS) residual values of the momentum and mass are set to be less than 10 -5 . The performance results are obtained after about 600 iterations.In the transient simulation, the boundary conditions is the same with that of the steady simulation. The interface between the rotor and stator is changed to "transient rotor stator". The accuracy of discretization in space is set as second order accuracy and second-order backward Euler scheme is used for the time discretization. The total calculating time for the pump is 12 revolutions and the every step time is 1.149×10-4s, corresponding to the rotating angle of 1°. Therefore, the results of final 3 revolutions are used to analyze the pressure fluctuation distributions in the pump. The steady simulation is used to calculate the performance of pump and the transient simulation is conducted to compare with the inner flow characteristics of the original pump and optimized one. The numerical simulations were carried out on a Dell Workstation with an Intel Xeon CPU E5-1650 (12 CPUs) and a clock speed of 3.2 GHz.

Fig. 6

 6 Fig.6 Flow chart of the optimization procedure

1 .

 1 Motor 2.Pump 3.Pressure sensor 4.Pressure sensor 5.Electromagnetic flowmeter 6.Valve 7.Valve 8.Tank

Fig 4 .Fig 5 .

 45 Fig 4. Diagram of the open test rig

7 .

 7 (a) Error analysis of approximation model of efficiency (b) Error analysis of approximation model of head Fig Error analysis of R square

  (a) Original scheme (b) Optimized scheme Fig 8. Comparison of velocity distribution in the impeller under design condition 3.4 Comparison of unsteady simulations In order to assess the unsteady pressure fluctuation distribution intensity directly and comprehensively, a statistical method is applied to analyze the pressure fluctuation intensity on each grid node in the fluid domains. Pressure standard deviation is calculated to indicate the pressure fluctuation intensity

  (a) Origin impeller (b) Optimized impeller Fig 9. Pressure fluctuation intensity distribution in the impeller

  (a) Orginal scheme (b) Optimized scheme Fig 10. pressure fluctuation intensity distribution in the volute (a) Time domain of pressure fluctuation (b) Frequency domain of pressure fluctuation Fig 11. Characteristics of pressure fluctuation of the monitor in the volute tongue

Table 1 .

 1 Main geometrical parameters of pump

		Geometrical parameters	Symbols	value
		Inlet diameter	Dj	50mm
		Outlet diameter	D2	160mm
		Blade outlet width	b2	6mm
		Blade inlet angel	β1	24°
	impeller	Blade outlet angel	β2	30°
		Blade wrap angle	φ	150°
		cutwater diameter	D3	170mm
		Inlet width	b3	18mm
	volute	Cutwater angle	α3	27°
		Outlet diameter	D4	32mm

Table 2

 2 Ranges of design parameters

	Design variables	Lower	Upper	origin
	b2/mm	5	8	6
	β2/°	27	33	30
	φ/°	135	165	150

Table 3

 3 Designed schemes and objective values

	Design schemes	b2 /mm	β2/°	φ/°	ηh/%	H/m
	1	7.84 28.26 146.05 71.40 9.20
	2	6.89 33.00 139.74 70.78 9.02
	3	6.42 29.84 157.11 75.27 8.82
	4	6.11 32.37 150.79 75.30 8.71
	5	7.37 32.68 152.37 71.93 9.30
	6	5.16 28.89 160.26 79.57 8.59
	7	5.47 29.53 147.63 78.20 8.71
	8	7. 21 31.11 165.00 73.12 9.21
	9	5.00 31.42 155.53 79.69 8.56
	10	5.79 29.21 135.00 76.95 8.77
	11	7.68 30.16 153.95 71.81 9.23
	12	7.05 28.58 136.58 73.47 9.10
	13	6.58 27.95 149.21 74.72 8.97
	14	5.32 27.00 142.89 78.93 8.63
	15	8.00 30.79 141.32 70.68 9.36
	16	7.53 27.63 158.68 72.82 9.22
	17	5.63 31.74 138.16 77.26 8.77
	18	5.95 32.05 163.42 75.96 8.80
	19	6.26 27.32 161.84 76.40 8.85
	20	6.74 30.47 144.47 72.64 8.89

Table 4

 4 Detailed parameters of multi-island GA

	Parameters of setup	values
	Interval of Migration	5
	Number of Generations	50
	Number of Islands	10
	Rate of Crossover	0.9
	Rate of Migration	0.01
	Rate of Mutation	0.01
	Sub-Population Size	10

Table 5

 5 Comparisons of performance between the origin and optimization

		b2/mm	β2/°	φ/°	η/% (Predicted)	H/m (Predicted)	η/%(CFD)	H/m(CFD)
	Original design	6	30	150	-	-	76.642	8.823
	Optimal design	5	29.64	156.1	79.89	8.54	79.845	8.838

ACKNOWLEDGMENTS

This study was supported by National Natural Science Foundation of China (Grant No. 51409123), Natural Science Foundation of Jiangsu Province Youth Fund (Grant No. BK20140554), China Postdoctoral Science Foundation (Grant No. 2014M560402) Postdoctoral Science Foundation of Jiangsu Province (Grant No. 1401069B) and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD). The author wishes to thank to the 4CPump research group.


Toru Shigemitsu

Hiroaki Fukuda

Junichiro Fukutomi

Wake and Potential Interference of Contra-Rotating Small-Sized Axial Fan at Design Flow Rate

Keywords: Small-sized axial fan, Contra-rotating rotors, Wake, Potential interference

Small-sized axial fans are used as air coolers for electric equipment. There is a strong demand for high pressure and large flow rate of fans according to the increase of quantity of heat from electric devices. Therefore, high rotational speed design is conducted, although, it causes the deterioration of the efficiency and increase of noise. Then, the adoption of contra-rotating rotors for small-sized fans is proposed for the improvement of the performance. In the case of the contra-rotating rotors, the wake and potential interference between the front and rear rotors influence on the performance and stable operation. Therefore, it is necessary to clarify the internal flow between the front and rear rotors. In the present paper, pressure fluctuations around the rotor phase locked each front and rear rotor's rotation are shown and the influences of the wake and potential interference at a design flow rate are discussed by the unsteady numerical analysis results.

INTRODUCTION

Because of spread of cloud computing, establishment of ubiquitous networking society and increase in rate of electric parts in machines, power consumption in data centers, IT devices and machines has been increasing significantly [1] . In the view of the issue of global warming and energy savings, there is a strong demand for the reduction of power consumption in above facilities and equipment. Electrical power used for the cooling of the IT devices for data centers is huge the same as that used for the IT devices itself in data centers. Small-sized axial fans are used as air coolers for electric equipment i.e. laptop, desk top computers and servers. There is a strong demand for high pressure and large flow rate of fans according to the increase of quantity of heat from electric devices. However, the increase of the pressure and flow rate by the increase of the fan diameter is restricted because of the limitation of the space. Therefore, high rotational speed design is conducted, although, it causes the deterioration of the efficiency and increase of the noise. On the other hand, lower rotational speed design [2] and advantages on the performance of the contra-rotating fans and pumps were verified by experimental results [START_REF] Furukawa | Experimental Study of Pump Characteristics of Contra-Rotating Axial Flow Pump[END_REF], [START_REF] Nouri | Design and Experimental Validation of a Ducted Counter-Rotating Axial-Flow Fans System[END_REF] . Then, the adoption of the contra-rotating rotors for the small-sized axial fans was proposed for the improvement of the performance. In the case of contra-rotating rotors, the axial space becomes larger than conventional small-sized axial fans. However, it is adequate choice to apply the contra-rotating rotors for small-sized fan because the axial space can be ensured in electrical devices as compared to that of the radial space.

In the case of the contra-rotating rotors, it is necessary to design the rear rotor considering the unsteady circumferential velocity distributions at the outlet of the front rotor [START_REF] Shigemitsu | Experimental Study on Rear Rotor Design in Contra-Rotating Axial Flow Pump[END_REF] . Further the fan noise becomes larger than the conventional rotor stator type fan because of the interaction between the front and rear rotors. Then, the passive noise reductions with the peforated blade was propsed for the contra-rotating fan [6] . It is important to clarify the influence of the wake from the front rotor to the rear rotor and the potential interaction between the front and rear rotors to increase the performance and to reduce the fan noise [START_REF] Sanders | Multi-Blade Row Interactions in a Transonic Axial Compressor: Part I -Stator Particle Image Velocimetry (PIV) Investigation[END_REF] . The blade row distance between front and rear rotors is a key parameter to consider the wake and potential interaction for the contra-rotating fan. The influence of the blade row distance between the front and rear rotors and pressure fluctuation on the casing wall were investigated for the counter rotating fan with fan diameter D=375mm [START_REF] Nouri | Design and Experimental Validation of a Ducted Counter-Rotating Axial-Flow Fans System[END_REF], [START_REF] Nouri | Experimental Study of the Instationary Flow Between Two Ducted Counter-Rotating Rotors[END_REF] . On the other hand, the conventional design method and the theory for the turbo machinery should be modified for small-sized axial fans because small-sized axial fans applied to electrical devices belong to extremely small size field in the turbo machinery [9] . Therefore, there is the strong demand to establish the design method for small-sized axial fans based on the internal flow between the front and rear rotors.

In the present paper, the performance curves of the contra-rotating small-sized axial fan with 100mm diameter are compared with the unsteady numerical analysis results to verify the validity of the unsteady numerical analysis results. After that, pressure fluctuations around the rotor phase locked each front and rear rotor's rotation are shown and the influences of the wake and the potential interference at the design flow rate are discussed by the unsteady numerical analysis results.

EXPERIMENTAL APPARATUS AND METHOD

The geometry of the rotor and the primary dimensions of the contra-rotating small-sized axial fan (RRtype) are shown in Fig. 1

N f =N r =1780[min -1 ]
. This fan was designed with the free vortex design. In this research, an aerofoil blade was used because there was a report that implied an advantage of the airfoil blade for the small-sized axial fan [10] , however a circular-arc blade was generally used for small-sized axial fans. Figure 2 shows the schematic diagram of the experimental apparatus for RRtype. The experimental apparatus was designed based on the Japanese Industrial Standard and the air blown in the test section passed the rotor, chamber, measurement duct and booster fan and blew out in the ambient atmosphere. The fan static pressure (P s ) was measured by the pressure difference between static holes downstream of the rotor installed at the chamber and ambient air. The fan static pressure of each front and rear rotor was also measured to investigate pressure curves of each front and rear rotor in wide flow rates range. Further, the rotational speed was controlled by the servo motor and the flow rates were measured by the orifice meter set at the measurement duct. The pressure curves of RRtype from the cutoff flow rate to the large flow rate were investigated in the experiment with the constant rotational speed N f =N r =1780[min -1 ]. In this research, the performance test was conducted with two different blade row distances L=10[mm] (non-dimensional axial distance divided by front rotor tip chord length L/l FR =0.42), 30[mm](L/l FR =1.27). L was defined as the axial distance from the trailing edge of the front rotor to the leading edge of the rear rotor at the hub.

NUMERICAL ANALYSIS CONDITIONS

The commercial software ANSYS-CFX 14.5 was used to investigate the flow condition and potential interference between the front and rear rotors of the contra-rotating small-sized axial fan which couldn't be measured by the experiment. In the numerical analysis, the numerical model which was almost the same with the experimental apparatus was used and three dimensional unsteady numerical analysis was conducted. The numerical analysis grids used for the numerical analysis are shown in Fig. 3. To simplify the numerical analysis, the servo motors and the shafts of the servo motors are removed in the numerical analysis model. The numerical domains comprised the inlet, rotor, chamber and outlet duct regions. The numerical grid elements were 617,090 for the inlet region, 1,479,336 for the chamber region and 237,628 for the outlet duct region. The numerical grid elements for the rotor region were 5,542,038 (L=10[mm]) or 6,038,826 (L=30[mm]) for RRtype. The tip clearance kept 1mm as the same with the experimental apparatus in the numerical analysis and the number of elements from the blade tip to the casing was 7. At the inlet boundary, the uniform velocity was given and the constant pressure was given at the outlet boundary condition. The coupling between the front and rear rotors was conducted by the sliding mesh (Transient Rotor Stator Model). The Near-Wall Treatment "Automatic" and SST (Shear Stress Transport) turbulence model were used. The unsteady numerical analysis was (corresponding 2° of the rotor rotation). The data of one rotor rotation were obtained after 6 rotor rotations in unsteady numerical analysis. The convergence criteria at each time step were set as the residual less than 1.0×10 -4 .

EXPERIMENTAL AND NUMERICAL ANALYSIS RESULTS

Performance curve of RRtype

The fan static pressure curves of RRtype at each blade row distance (L=10, 30[mm]) are shown in Fig. 4. In Fig. 4, the fan static pressure obtained by both experiment and numerical analysis are shown. In the unsteady numerical analysis, the fan static pressure was obtained at the same point in the experiment and 180 static pressure data of one rotor rotation were averaged. The vertical axis and horizontal axis in Fig. 4 show the fan static pressure P s and flow rate Q. It was found from Fig. 4 that fan static pressure of RRtype linearly increased as the flow rate decreased, and the pressure curve of RRtype showed the stable negative curve. The fan static pressure curves of the front rotor and rear rotor are shown in Fig. 5. The fan static pressure of the front rotor also showed the same tendency of RRtype, but the fan static pressure of the rear rotor represented a different tendency from the fan static pressure of RRtype in the partial flow rates. And the slight positive slope of the pressure curve was confirmed in the low flow rates for the rear rotor. It was clarified from the previous research [11] that the fan static pressure kept almost constant by the blade row distance L=30[mm] (L/l FR =1.27). It was confirmed from Fig. 4 that the difference of the fan static pressure between L=10[mm] and L=30[mm] was P s =0.02 (Exp.) at the design flow rate Q d =0.016[m [START_REF] Furukawa | Experimental Study of Pump Characteristics of Contra-Rotating Axial Flow Pump[END_REF] /s]. The influence of the blade row distance in the range of L=10-30[mm] on the fan static pressure was small. The numerical analysis results could predict the experimental results accurately and capture the tendency of the performance curves of the experimental results. It was difficult to measure the efficiency of this small fan. Therefore, the efficiency of the test fan was investigated by the numerical analysis results. The static pressure efficiency of RRtype for L=10[mm] and 30[mm] is shown in Fig. 6. The maximum static pressure efficiency was η=58.0% for L=10[mm] and η=57.0% for L=30[mm] and the difference of the efficiency for L=10[mm] and 30 [mm] were small in wide flow rate range. The static pressure efficiency of each front and rear rotors are shown in Fig. 7. The static pressure efficiency of the rear rotor was larger than that of the front rotor especially around the design flow rate. The static pressure efficiency of the front rotor for L=30[mm] was larger than that for L=10[mm] at the design flow rate and vice versa in the case of the rear rotor. It is important to clarify the wake and potential interference for the stable operation of the contra-rotating small sized fan, therefore, the interferences at the design flow rate were investigated by the unsteady numerical analysis.

Q[m 3 /s] Front Rear L=10[mm] L=30[mm]

Wake and potential interference between front and rear rotors (L=10, 30[mm])

In this section, the pressure fluctuations around the rotor phase locked each front and rear rotor's rotations are shown and the influences of the wake and potential interference are discussed by the unsteady numerical analysis results. The pressure fluctuations at the LE and the middle position of the SS of the rear rotor obtained by the numerical analysis are shown in Figs. 11 and12 respectively. The flow rate, vertical axis and horizontal axis are the same in Figs.9 and 10. Because the pressure in Figs. 11 and12 were phase locked the rear rotor, the wake and potential interference from the front rotor became clear if the unsteady pressure fluctuations from the rear rotor were small. The periodic pressure fluctuations with four peaks corresponding to the front rotor blade number Z F =4 were generated for the blade row distance L=10mm and L=30mm in Figs. 11 and12. The pressure fluctuations ∆P s at the LE of the rear rotor were larger than that at the TE of the front rotor. The pressure fluctuations ∆P s at the LE and middle position of the SS of the rear rotor near the shroud at r/r c =0.98 were significantly large compared to the other radial positions. These large pressure fluctuations near the shroud would be related to the flow condition at the rear rotor inlet. When the pressure fluctuations ∆P s of the front and rear rotors were compared, the pressure fluctuations ∆P s of the rear rotor larger than that of front rotor as shown in Figs.9-12. The variation of the blade loading of the rear rotor with the rotor rotation caused by the pressure fluctuations in Figs.11 and 12 would be larger than that of the front rotor. Therefore, it is important to consider the influence of the wake and potential interferences especially on the rear rotor for the stable operation of the contra-rotating small sized fan. On the other hand, the wake and potential interference from the front rotor could be suppressed with the increase of the blade row distance from L=10mm to L=30[mm] as shown in Figs. 11 and12.

The pressure fluctuations of the rear rotor were large. Then, the attack angle variations with the rotor rotation at the design flow rate Q d were obtained by the numerical analysis. Figure 13 shows the attack angle variations of the rear rotor with the rotors rotation. The vertical axis is attack angle and the horizontal axis is the relative angle of the front and rear rotors θ fr . The attack angle variation of the rear rotor for L=10 [mm] were larger than that for L=30 [mm]. It was considered that the wake from the front rotor would influence on the attack angle of the rear rotor. The variation of the attack angle might induce the pressure fluctuations in θ rf [deg] (a) L=10mm appropriate blade row distance because the performance decrease with the increase of the blade row distance was small by L=30mm and the wake and potential interference could be suppressed effectively for L=30mm.

CONCLUDING REMARKS

The performance test of the contra-rotating small-sized axial fan was conducted under the condition of different blade row distance. Then, the wake and potential interference between the front and rear rotors were investigated at the design flow rate Q=0.016[m [START_REF] Furukawa | Experimental Study of Pump Characteristics of Contra-Rotating Axial Flow Pump[END_REF] /s] with the unsteady numerical analysis. As a result, following concluding remarks could be obtained. 3. The attack angle variation of the rear rotor for L=10 [mm] were larger than that for L=30 [mm]. It was considered that the wake from the front rotor would influence on the attack angle of the rear rotor. The blade row distance L=30[mm] would be appropriate blade row distance because the performance decrease with the increase of the blade row distance was small by L=30[mm] and the wake and potential interference could be suppressed effectively for L=30[mm].

Figure. 2 Figure. 1 Figure. 3

 213 Figure. 2 Experimental apparatus

Figure. 4 Figure. 5 Figure. 6 Figure. 7

 4567 Figure. 4 Fan static pressure of RRtype

  The sampling points for pressure fluctuations and definition of the relative position of the front and rear rotors are shown in Fig.8. The respective data for the front rotor were obtained at the middle position of the pressure surface (PS) of the blade and 0.1[mm] downstream from the trailing edge (TE) of the front rotor. And the data for the rear rotor were obtained at the middle position of the suction surface (SS) of the blade and 0.1mm upstream from the leading edge (LE) of the rear rotor. The rotational angles of the front and rear rotors are shown as θ f , θ r respectively. The rotational angles were defined as the angles of the 25% chord length position of each front and rear rotor from the arbitrary meridional plane. The relative angle of the front and rear rotors θ fr are defined as θ fr =θ f +θ r , the initial relative angle is θ fr =0[deg.] (θ f =θ r =0[deg.]). The rotational direction of the angle θ f, θ r were the same with the rotational direction of the front and rear rotors respectively. The pressure fluctuations at the middle position of the PS and the TE of the front rotor at the design flow rate Q d =0.016[m 3 /s] obtained by the numerical analysis are shown in Figs.9 and 10 respectively. The vertical axis is the fluctuation of the static pressure ∆P s from the averaged value and the horizontal axis is the relative angle of the front and rear rotors θ fr . Further, r/r c shows non-dimensional radius divided by the radius at the casing; r/r c =0.45 and r/r c =1.0 correspond the hub and casing. Because the pressure in Figs.9 and 10 were phase locked the front rotor, the potential interference from the rear rotor became clear if the unsteady pressure fluctuations from the front rotor were small. The periodic fluctuations with five peaks corresponding to the rear rotor blade number Z R =5 were generated for the blade row distance L=10mm and L=30mm in Figs.9 and 10. The pressure fluctuations ∆P s at the middle position of the PS of the front rotor for L=10mm and L=30mm were within ± 0.75Pa (fan static pressure P s =11.47Pa at Q d for L=10mm) and ±0.25Pa (fan static pressure P s =11.45Pa at Q d for L=30[mm]) respectively. On the other hand, the pressure fluctuations ∆P s at the TE of the front rotor for L=10mm and L=30[mm] were within ±1.2[Pa] and ±0.3[Pa] respectively. The pressure fluctuation became small as the distance from the rear rotor became large and the potential interference of the rear rotor could be suppressed with the increase of the blade row distance from L=10[mm] to L=30[mm].

Figure. 8 1 -Figure. 9

 819 Figure. 8 Sampling points and relative position of front and rear rotors

Fig. 11 .Figure. 10

 1110 Figure. 10 Pressure fluctuations at TE of front rotor

Figure. 11 Figure. 12

 1112 Figure. 11 Pressure fluctuations at LE of rear rotor (b) L=30mm

1 .

 1 The difference of the fan static pressure between L=10[mm] and L=30[mm] was P s =0.02 (Exp.) at the design flow rate Q d =0.016[m 3 /s]. The maximum static pressure efficiency was η=58.0% for L=10[mm] and η=57.0% for L=30[mm]. The influences of the blade row distance in the range of L=10-30[mm] on the fan static pressure and static pressure efficiency were small. 2. The periodic pressure fluctuations corresponding to the blade number of each front and rear rotor were confirmed at the LE of the rear rotor and TE of the front rotor respectively. The potential interference from the rear rotor could be suppressed with the increase of the blade row distance from L=10[mm] to L=30[mm]. The wake and potential interference from the front rotor could be suppressed also with the increase of the blade row distance from L=10[mm] to L=30[mm].

Figure. 13

 13 Figure. 13 Attack angle variations of rear rotor with rotation (b) L=30mm

  and Table 1, respectively. The hub tip ratio was D h /D t =45[mm]/98[mm], tip clearance was c=1[mm] and the design flow rate was Q d =0.016[m 3 /s]. Fan static pressure at the design point was P dRR =14.7[Pa] for RRtype with the same fan static pressure of each front and rear rotor. The rotational speed of front and rear rotors of RRtype was

Table . 1

 . Primary dimensions of RRtype

			Rear rotor
				Front rotor
			Hub	Mid Tip
		Diameter[mm]	45	72	98
		Blade Number		4
	Front Rotor	Blade Profile Solidity	NACA 4409 1.196 0.496 0.29
		Stagger Angle	44.7° 61.1° 68.2°
		Blade Number		5
	Rear Rotor	Blade Profile Solidity	NACA 4412 0.91 0.447 0.288
		Stagger Angle	56.7° 64.5° 69.6°

ACKNOWLEDGMENTS

The authors wish to show our special thanks to the supports by the Komiya research aid, the project research aid from The University of Tokushima and Japan Science and Technology Agency.

for Gas Turbines and Power, Vol.135, No.2, 022601, doi:10.1115/1.4007756, 2013.


Debayan Maity
email: debayan.maity@u-bordeaux.fr

AND Marius Tucsnak
email: marius.tucsnak@u-bordeaux.fr

Enrique Zuazua
email: enrique.zuazua@deusto.es

CONTROLLABILITY AND POSITIVITY CONSTRAINTS IN POPULATION DYNAMICS WITH AGE STRUCTURING AND DIFFUSION

Keywords: Population dynamics, Null controllability. AMS subject classifications. 93B03, 93B05, 92D25

In this article, we study the null controllability of a linear system coming from a population dynamics model with age structuring and spatial diffusion (of Lotka-McKendrick type). The control is localized in the space variable as well as with respect to the age. The first novelty we bring in is that the age interval in which the control needs to be active can be arbitrarily small and does not need to contain a neighbourhood of 0. The second one is that we prove that the whole population can be steered into zero in a uniform time, without, as in the existing literature, excluding some interval of low ages. Moreover, we improve the existing estimates of the controllability time and we show that our estimates are sharp, at least when the control is active for very low ages. Finally, we show that the system can be steered between two positive steady states by controls preserving the positivity of the state trajectory. The method of proof, combining final-state observability estimates with the use of characteristics and with L ∞ estimates of the associated semigroup, avoids the explicit use of parabolic Carleman estimates.

Introduction and main results

In this article, we study the null-controllability of an infinite dimensional linear system describing the dynamics of a single species age-structured population with spatial diffusion. In these models, going back to Gurtin [START_REF] Gurtin | A system of equations for age-dependent population diffusion[END_REF] and generalizing the classical Lotka-McKendrick system, the state space of the system is H = L 2 ([0, a † ] × Ω), where a † denotes the maximal age an individual can attain and Ω ⊂ R n (with n ∈ N in general but with n = 3 for real life applications) is an open bounded set which represents the spatial environment occupied by the individuals. Let p(t, a, x) be the distribution density of individuals with respect to age a 0 and spatial position x ∈ Ω at some time t 0. Then, according to the above reference, the function p satisfies the degenerate parabolic partial differential equation ∂p ∂t + ∂p ∂a -Lp + µ(a)p = mv (t, a, x) ∈ (0, ∞) × (0, a † ) × Ω, (

where the operator L is defined by

Lp = n i,j=1 ∂ ∂x i σ ij ∂p ∂x j , (1.2)
Date: October 4, 2018. Debayan Maity and Marius Tucsnak acknowledge the support of the Agence Nationale de la Recherche -Deutsche Forschungsgemeinschaft (ANR -DFG), project INFIDHEM, ID ANR-16-CE92-0028. The research of Enrique Zuazua was supported by the Advanced Grant DyCon (Dynamical Control) of the European Research Council Executive Agency (ERC), the MTM2014-52347 and MTM2017-92996 Grants of the MINECO (Spain) and the ICON project of the French ANR-16-ACHN-0014.

with σ ij = σ ji ∈ C 2 (Ω), for 1 i, j n, and we assume there exists a constant c > 0 such that n i,j=1

σ ij (x)ξ i ξ j c|ξ| 2 (x ∈ Ω, ξ ∈ R n ).
Moreover, the positive function µ denotes the natural mortality rate of individuals of age a, supposed to be independent of the spatial position x and of time. The control function is v, depending on t, a and x, whereas m is the characteristic function of (a 1 , a 2 ) × ω, with 0 a 1 < a 2 a † and ω ⊂ Ω an open set. Thus the control is localized both in age and with respect to the spatial variable. This control process corresponds to harvesting of adding individuals of age between a 1 and a 2 from the spatial domain ω. Note that equation (1.1) is a slight generalization of the one proposed in [START_REF] Gurtin | A system of equations for age-dependent population diffusion[END_REF], where the operator L is just the standard Laplacian. We denote by β the positive function describing the fertility rate at age a, supposed to be independent of the spatial position x and of time, so that the density of newly born individuals at the point x at time t is given by p(t, 0, x) = a † 0 β(a)p(t, a, x) da (t, x) ∈ (0, ∞) × Ω. (1.3) We assume that the individuals never leave the set Ω, so that p satisfies the Neumann boundary condition

∂p ∂ν L = n i,j=1
σ ij ∂p ∂x j n i = 0 (t, a, x) ∈ (0, ∞) × (0, a † ) × ∂Ω, (1.4) where n denotes the unit outer normal to ∂Ω. To complete the model, we introduce the initial condition p(0, a, x) = p 0 (a, x) (a, x) ∈ (0, a † ) × Ω.

(

We assume that the fertility rate β and the mortality rate µ satisfy the conditions (H1) β ∈ L ∞ (0, a † ), β 0 for almost every a ∈ (0, a † ). (H2) µ ∈ L 1 [0, a * ] for every a * ∈ (0, a † ), µ 0 for almost every a ∈ (0, a † ).

(H3) a † 0 µ(a) da = +∞.

For more details about the modelling of such system and the biological significance of the hypotheses, we refer to Webb [START_REF] Webb | Theory of nonlinear age-dependent population dynamics[END_REF].

Theorem 1.1. Assume that β and µ satisfy the conditions (H1)-(H3) above. Moreover, suppose that the fertility rate β is such that β(a) = 0 for all a ∈ (0, a b ), (1.6) for some a b ∈ (0, a † ) and that a 1 < a b . Recall that m is the characteristic function of (a 1 , a 2 ) × ω, with 0 a 1 < a 2 a † and that ω ⊂ Ω is an open set. Then for every τ > a 1 + a †a 2 and for every

p 0 ∈ L 2 ((0, a † ) × Ω) there exists a control v ∈ L 2 ((0, τ ) × (a 1 , a 2 ) × ω) such that the solution p of (1.1)-(1.5) satisfies p(τ, a, x) = 0 for all a ∈ (0, a † ), x ∈ Ω. (1.7)
To state our result on controllability with positivity constraints, we first define the concept of non-negative steady state for (1.1) - (1.5).

Definition 1.2. Let v s ∈ L ∞ ((0, a † ) × Ω) be a steady interior control such that v s 0 a.e. on (0, a † ) × Ω. A non-negative function p s ∈ L ∞ ((0, a † ) × Ω) satisfying the equations                ∂p s ∂a -Lp s + µ(a)p s = mv s (a, x) ∈ (0, a † ) × Ω, ∂p s ∂ν L = 0 (a, x) ∈ (0, a † ) × ∂Ω, p s (0, x) = a † 0 β(a)p s (a, x) da, x ∈ Ω, (1.8)
is said to be a non-negative steady state for (1.1) -(1.5).

Our second main result can be stated as follows:

Theorem 1.3. Assume the hypothesis of Theorem 1.1. Let p s,I and p s,F are two non-negative steady states of the system (1.1) -(1.5). Assume that there exist a * ∈ (0, a † ) and δ > 0 such that p s,I (a, x), p s,F (a, x) δ a.e. on [0, a * ] × Ω.

(1.9)

Then there exist τ > 0 and v ∈ L ∞ ((0, τ ) × (0, a † ) × Ω) such that the problem (1.1) -(1.5) with p 0 (a, x) = p s,I (a, x)
admits a unique solution p satisfying p(τ, a, x) = p s,F (a, x) for all (a, x) ∈ (0, a † ) × Ω.

Moreover, p(τ, a, x) 0 for a.e. (t, a, x) ∈ (0, τ ) × (0, a † ) × Ω. • if R < 1 then there exists a unique non-negative solution to (1.8)

• if R = 1 and v s ≡ 0, then there exists infinitely many solutions to (1.8) of the form p s = αe -a 0 µ(r) dr , α ∈ (0, ∞). • if R > 1 then there is no non-negative solution to (1.8).

Consequently, the existence of non-negative steady states satisfying (1.9) is ensured at least when R = 1 and v s = 0. Another situation where we know that such states exist is when the control is active on all [0, a † ] × Ω and R < 1 (see [START_REF] Ainseba | Local exact controllability of the age-dependent population dynamics with diffusion[END_REF]Theorem 3.1]).

On the other hand, one can show that, if p s is a non-negative solution to the system (1.8), then

lim a→a † p s (a, x) = 0 a.e. x ∈ Ω. (1.10)
Therefore we cannot assume that the initial or target non-negative steady states are bounded from below by a strictly positive constant on [0, a † ], i.e., we cannot take a * = a † in Theorem 1.3.

Let us now mention some related works from the literature. The null controllability results of the diffusion free age-dependent population dynamics model were first obtained by Barbu, Iannelli and Martcheva [START_REF] Barbu | On the controllability of the Lotka-McKendrick model of population dynamics[END_REF]. They proved the state of the system can be steered to any steady state, except for a small interval of ages near zero. Recently, Hegoburu, Magal and Tucsnak [START_REF] Hegoburu | Controllability with positivity constraints of the lotka-mckendrick system[END_REF] proved that this restriction is not necessary, provided individuals do not reproduce at the age close to zero. They also proved there exists controls which preserves the positivity of the state trajectory. However, in the both works, the control is supported in the interval (0, a 0 ), for some a 0 < a † . Recently, Maity [START_REF] Maity | On the Null Controllability of the Lotka-Mckendrick System[END_REF] proved that null controllability can be achieved by controls supported in any subinterval of (0, a † ), provided we control before the individuals start to reproduce.

Concerning the models with spatial diffusion, namely for the system (1.1) -(1.5), as far as we know, the first result was obtained by Ainseba and Anit ¸a [START_REF] Ainseba | Local exact controllability of the age-dependent population dynamics with diffusion[END_REF]. They proved that the system (1.1) -(1.5) can be driven to a steady state in any arbitrary time τ > 0 keeping the positivity of the trajectory, provided the initial data is close to the steady state and the control acts in a spatial subdomain ω ⊂ Ω but for all ages. When control acts in a spatial subdomain and only for small ages, a similar result for a large time was proved by Ainseba and Anit ¸a [START_REF] Ainseba | Internal exact controllability of the linear population dynamics with diffusion[END_REF]. In [START_REF] Ainseba | Exact and approximate controllability of the age and space population dynamics structured model[END_REF] Ainseba proved null controllability of the system (1.1)-(1.5) except for a small interval of ages near zero, with controls acting everywhere in the ages but in a spatial subdomain. Recently, Hegoburu and Tucsnak [START_REF] Hegoburu | Null controllability of the Lotka-McKendrick system with spatial diffusion[END_REF] proved that the system (1.1) -(1.5) is null controllable for all ages and in any time by controls localized with respect to the spatial variable but active for all ages. Their method is based on Lebeau-Robbiano type strategy, originally developed for the null-controllability of the heat equation. Traore [START_REF] Traore | Null controllability of a nonlinear population dynamics problem[END_REF] considered a similar model with nonlinear distributions of the newborns. He proved null controllability except for small ages with controls localized in space variable and active for all ages. Martinez et. al [START_REF] Martinez | Regional null controllability of a linearized Crocco-type equation[END_REF] considered linearized Croco-type equation, which is similar to the system (1.1) -(1.5), with β = µ = 0. They proved regional null controllability of such system.

The main novelties brought in by our paper are:

• We improve the existing estimates on the time necessary to control the system to zero and we show that our global controllability result applies to individuals of all ages, without needing to exclude ages in a neighbourhood of zero. • We are able to tackle the case of a control which is active for ages a ∈ [a 1 , a 2 ], with arbitrary a 1 ∈ [0, a † ) and a 2 ∈ (a 1 , a † ], provided that supp β ∩ [0, a 1 ] = ∅. Thus, unlike in the existing literature, we do not need to control arbitrarily low ages. • Unlike most of the approaches in the literature, our methodology does not require adaptations of the existing parabolic Carleman estimates to the adjoint system of (1.1) -(1.5). We just combine characteristics method with existing observability estimates for parabolic equations. Thus our approach applies independently of the method used to derive final-state observability for the associated parabolic system (moment methods, local or global Carleman estimates, Lebeau-Robbiano strategy,...). • Controllability with positivity constraints is proved, as far as we know for the first time, with a control which is localized both in age and with respect to the space variable. The methodology employed to obtain this result is based on duality and L ∞ estimates for parabolic PDEs. The remaining part of this work is organized as follows:

• In Section 2 we first recall some basic facts about the Lotka-McKendrick semigroup with diffusion. We next formulate our control problem in a semigroup setting and we define the associated adjoint semigroup. • In Section 3 we prove the final state observability for the adjoint system and, as a consequence, we obtain the proof of the main result in Theorem 3.2. • Section 4 is devoted to the proof that controllability between positive steady states can be achieved in sufficiently large time, i.e., to the proof of Theorem 1.3. • Section 5 is devoted to the description of possible extensions and open questions.

Lotka-McKendrick Semigroup with Diffusion

In this section, we provide some basic results on the population semigroup for the linear age structured model with diffusion and its adjoint operator. Most of them were existing in the literature, so we just give the statements and the appropriate references. In some cases, namely when the adjoint operator is involved, we did not find detailed justifications in the existing literature, so, with no claim of originality, we felt necessary to give a more detailed presentation.

We write below equations (1.1)-(1.5) as an abstract control system with input space

H = L 2 (0, a † ; L 2 (Ω))
Before introducing the semigroup generator, we consider the diffusion free population operator

A 1 : D(A 1 ) → H, defined by D(A 1 ) = ϕ ∈ H | ϕ(•, x) is locally absolutely continuous on [0, a † ), ϕ(0, x) = a † 0 β(a)ϕ(a, x) da for a.e. x ∈ Ω, ∂ϕ ∂a + µϕ ∈ H , A 1 ϕ = - ∂ϕ ∂a -µϕ, (2.1)
and the diffusion operator

A 2 : D(A 2 ) → H defined by D(A 2 ) = ϕ ∈ L 2 (0, a † ; H 2 (Ω)) | ∂ϕ ∂ν L = 0 on ∂Ω , A 2 ϕ = n i,j=1 ∂ ∂x i σ ij ∂ϕ ∂x j (ϕ ∈ D(A 2 )). (2.2)
We also introduce the input space U = H and the control operator B ∈ L(U, H) defined by

Bu = mu (u ∈ U ). (2.3)
With the above notation, we rewrite the system (1.1)-(1.5) as:

ż(t) = Az(t) + Bu(t), (2.4)
z(0) = p 0 , (2.5
) where we have set p(t, •) = z(t), v(t, •) = u(t) and the population operator with diffusion A : D(A) → H is defined by

D(A) = D(A 1 ) ∩ D(A 2 ), A = A 1 + A 2 .
(2.6) The fact that the system we consider is well-posed follows from the following result: Lemma 2.1. The operator (A, D(A)) is the infinitesimal generator of a strongly continuous semigroup T on H. Proof. A proof of this lemma can be found in [START_REF] Huyer | Semigroup formulation and approximation of a linear age-dependent population problem with spatial diffusion[END_REF]Theorem 2.8].

With the above notation, our main result in Theorem 1.1 can be rephrased to : the pair (A, B), with A defined in (2.6) and B defined in (2.3) is null controllable in ant time τ > a 1 +max{a 1 , a † -a 2 }. It is well-known that, the null controllability in time τ of (A, B) is equivalent to the final-state observability in time τ of the pair (A * , B * ), where A * and B * are the adjoint operators of A and B, respectively (see, for instance, [START_REF] Tucsnak | Observation and control for operator semigroups[END_REF]Section 11.2]). It is thus important to determine the adjoint of the operator A.

To this aim, we introduce an auxiliary unbounded operator (A 0 , D(A 0 )) defined by

D(A 0 ) = ϕ ∈ H | ϕ(•, x) is locally absolutely continuous on [0, a † ), ϕ ∈ L 2 (0, a † ; H 2 (Ω)), lim a→a † ϕ(a, x) = 0 for a.e. x ∈ Ω, ∂ϕ ∂ν L = 0 on ∂Ω, ∂ϕ ∂a -µϕ + Lϕ ∈ H , A 0 ϕ = ∂ϕ ∂a -µϕ + Lϕ (ϕ ∈ D(A 0 )). (2.7)
We have the following lemma:

Lemma 2.2. The operator (A 0 , D(A 0 )) is the infinitesimal generator of a strongly continuous semigroup T 0 on H.

Proof. For ϕ ∈ D(A 0 ), we have

(A 0 ϕ, ϕ) H = lim a→a - † a 0 Ω ∂ϕ ∂a -µϕ + Bϕ ϕ = lim a→a - † Ω ϕ 2 (a) 2 - Ω ϕ 2 (0) 2 -lim a→a - † a 0 Ω µϕ 2 + a ij ∂ϕ ∂x i ∂ϕ ∂x j 0, thus A 0 is a dissipative operator.
To show that it is m-dissipative, we prove below that I -A 0 is onto. To this aim, let f ∈ H and consider the equation, of unknown ϕ ∈ D(A 0 ),

ϕ - ∂ϕ ∂a + µϕ -Lϕ = f in (0, a † ) × Ω, lim a→a - † ϕ(a, x) = 0, ∂ϕ ∂ν L = 0 on ∂Ω. (2.8) Denoting ϕ(a, x) = exp -a - a 0 µ(r) dr ϕ(a, x) (a ∈ (0, a † ), x ∈ Ω), (2.9)
we see that (2.8) is equivalent to the equation, of unknown ϕ,

- ∂ ϕ ∂a -L ϕ = f in (0, a † ) × Ω, ϕ(a † , x) = 0 in Ω, ∂ ϕ ∂ν L = 0 on ∂Ω, where f (a, x) = exp -a - a 0 µ(r) dr f (a, x). It is easy to see that f ∈ H. It is easily seen that the above equation has a solution ϕ ∈ L 2 (0, a † ; H 2 (Ω)) ∩ H 1 (0, a † ; L 2 (Ω)) with ϕ(a, •) L 2 (Ω) C -a - a 0 µ(r) dr f H .
The above estimate and (2.9) imply that (2.8) has a solution ϕ ∈ D(A 0 ), thus A 0 is m-disiipative.

Hence A 0 generates a C 0 semigroup on H.

We are now in a position to rigourously construct the adjoint of the unbounded operator A.

Proposition 2.3. The adjoint of (A, D(A)) in H is defined by

D(A * ) = D(A 0 ), A * ψ = ∂ψ ∂a -µψ + βψ(0, x) + Lψ. (2.10)
Proof. We can easily verify that D(A 0 ) ⊂ D(A * ).

To prove the reverse inclusion, we first note that for λ > 0 large enough, the operator λI -A 0 is boundedly invertible (this follows from the that A 0 is a semigroup generator). For those values of λ we can thus consider the operator I -F(λ), where F(λ) ∈ L(L 2 (Ω)) is defined by

F(λ)g(x) = (λI -A 0 ) -1 (β(a)g(x)) (0, x).
(2.11)

We note that for λ large enough the operator I -F(λ) is invertible. Indeed, this follows from the fact that lim λ→∞

F(λ) L(L 2 (Ω),H) = 0.
For λ as above, we define G λ : H → H defined by by G λ f = ϕ λ where ϕ λ solves

λϕ λ - ∂ϕ λ ∂a + µϕ λ -Lϕ λ -β(a)ϕ λ (0, x) = f in (0, a † ) × Ω, lim a→a - † ϕ λ (a, x) = 0, ∂ϕ λ ∂ν L = 0 on ∂Ω. (2.
12)

The fact that the operator G λ is well defined follows from the fact that the unique solution of (2.12) is clearly given by ϕ λ (a, x) = (λI

-A 0 ) -1 (f (a, x) + V λ,f (a, x)) , (2.13) where V λ,f (a, x) = β(a) (I -F (λ)) -1 (λI -A 0 ) -1 f (0, x) .
This means, in particular, that ϕ λ ∈ D(A 0 ). We are now in the position to prove the inclusion D(A * ) ⊂ D(A 0 ). To this aim, take λ as above and let ψ ∈ D((λI -A) * ). Then there exists f ∈ H such that

a † 0 Ω ψ(λI -A)ϕ = a † 0 Ω
f ϕ for all ϕ ∈ D(A).

Let η λ = G λ f , with G λ defined several lines above. Then, using (2.12) and integrating by parts we obtain

a † 0 Ω f ϕ = lim a→a - † a 0 Ω λη λ - ∂η λ ∂a + µη λ -Lη λ -β(a)η λ (0, x) ϕ = lim a→a - † Ω a 0 η λ λϕ + ∂ϕ ∂a + µϕ -Lϕ = a † 0 Ω η λ (λI -A)ϕ.
Therefore,

a † 0 Ω (ψ -η λ )(λI -A)ϕ = 0, for all ϕ ∈ D(A). (2.14)
By choosing ϕ = (λI -A) -1 (ψ -η λ ) we get a † 0 Ω |ψ -η λ | 2 = 0.
Thus ψ ∈ D(A 0 ) and ψ solves (2.12). This completes the proof of the proposition.

An Observability Inequality:

As mentioned above, the null-controllability of a pair (A, B) is equivalent to the final state observability of the pair (A * , B * ), see [START_REF] Tucsnak | Observation and control for operator semigroups[END_REF]Theorem 11.2.1]. Recall that that final-state observability of (A * , B * ) is defined as

Definition 3.1. [19, Definition 6.1.1] The pair (A * , B * ) is final state observable in time τ if there exists a k τ > 0 such that τ 0 B * T * t q 0 H k 2 τ T * τ q 0 2 (q 0 ∈ D(A * )).
For A defined in (2.6) and q 0 ∈ H we set

q(t) = T * t q 0 (t 0),
where T is the semigroup generated by A. According to Proposition 2.3 we have:

                 ∂q ∂t - ∂q ∂a -Lq -β(a)q(t, 0, x) + µ(a)q = 0, t 0, (a, x) ∈ (0, a † ) × Ω, q(t, a † , x) = 0, t 0, x ∈ Ω, ∂q ∂ν L = 0, t 0, (a, x) ∈ (0, a † ) × ∂Ω, q(0, a, x) = q 0 (a, x), (a, x) ∈ (0, a † ) × Ω. (3.1)
In view of [START_REF] Tucsnak | Observation and control for operator semigroups[END_REF]Theorem 11.2.1], the statement in Theorem 1.1 is equivalent to the following theorem:

Theorem 3.2. Under the assumption of Theorem 1.1, the pair (A * , B * ) is final-state observable for every τ > a 1 + a †a 2 . In other words, for every τ > a 1 + a †a 2 there exists k τ > 0 such that the solution q of (3.1) satisfies

a † 0 Ω q 2 (τ, a, x) dxda k 2 τ τ 0 a 2 a 1 ω q 2 (t, a, x) dxdadt (q 0 ∈ D(A * )). (3.2)
Before we begin the proof of the above result, let us briefly describe its main steps. The first one is writing an explicit expression of q. We define

V (t, x) = q(t, 0, x), (t, x) ∈ (0, τ ) × Ω. (3.3)
Integrating along the characteristic lines, the solution of (3.1) can be written as

q(t) =          π(a) π(a + t) e tA 2 q 0 (a + t, •) + t 0 π(a) π(a + t -s) e (t-s)A 2 β(a + t -s)V (s, •) t a † -a, t t+a-a † π(a) π(a + t -s) e (t-s)A 2 β(a + t -s)V (s, •) ds t > a † -a, (3.4)
where π(a) = e -a 0 µ(r)dr . Without loss of generality, we can assume that a 2 a b and let τ be as in Theorem 3.2. We decompose the interval (0, a † ) as

(0, a † ) = (0, a) ∪ ( a, a † ),
where a 0 is chosen suitably so that a 0 < a 2 and τ > a † -a for all a ∈ (a 0 , a † ). Now using the expression of q in (3.4) and choosing a 0 suitably, we can show that

a † 0 Ω q 2 (τ, a, x) dxda = a 0 0 Ω q 2 (τ, a, x) dxda + a † a 0 Ω q 2 (τ, a, x) dxda (3.5) C a 0 0 Ω q 2 (τ, a, x) dxda + τ η Ω q 2 (t, 0, x) dxdt , (3.6)
for some η > a 1 (see proof of Theorem 3.2 for more details). The second step consists in deriving upper appropriate bounds for each one of the terms in the right-hand side of (3.5). This is accomplished by combining some change of variables using the characteristics of the diffusion free problem with some known observability inequalities for parabolic equations.

To accomplish this program, we first recall the following observability inequality for parabolic equations (see, for instance, Imanuvilov and Fursikov [START_REF] Fursikov | Controllability of evolution equations[END_REF]) : Proposition 3.3. Let T > 0, 0 t 0 < τ and t 1 ∈ (t 0 , T ]. Then for every w 0 ∈ L 2 (Ω), the solution w of the initial and boundary problem

           ∂w ∂s (s, x) -Lw(s, x) = 0 ((s, x) ∈ (t 0 , T ) × Ω), ∂w ∂ν L = 0 ((s, x) ∈ (t 0 , T ) × ∂Ω), w(t 0 , x) = w 0 (x), (x ∈ Ω), (3.7)
satisfies the estimate

Ω w 2 (T, x) dx Ω w 2 (t 1 , x) dx c 1 e c 2 t 1 -t 0 t 1 t 0 ω w 2 (s, x) dx ds, (3.8)
where the constants c 1 and c 2 depend on L, on Ω and on τ .

In the following two propositions we estimate each of the terms appearing in the right-hand side of (3.5).

Proposition 3.4. Let us assume the hypothesis of Theorem 1.1 and let τ > a 1 and a 0 ∈ (0, a b ). Then there exists a constant C > 0 such that, for every q 0 ∈ D(A * ), the solution q of the system (3.1), obeys

a 0 0 Ω q 2 (τ, a, x) dxda C τ τ 0 a 2 a 1 ω q 2 (t, a, x) dx da dt. (3.9)
Proof. First of all, without loss of of generality we can assume that a 2 a b (otherwise we simply observe for small ages). We can also assume that a 0 > max{a 1 , a 2τ }. Since β(a) = 0 for all a ∈ (0, a b ), q satisfies

       ∂q ∂t - ∂q ∂a -Lq + µ(a)q = 0, t 0, (a, x) ∈ (0, a b ) × Ω, ∂q ∂ν L = 0 t 0, (a, x) ∈ (0, a b ) × ∂Ω. (3.10)
This means that q satisfies the adjoint of a system called "Crocco type" (see [START_REF] Martinez | Regional null controllability of a linearized Crocco-type equation[END_REF]), where the authors proved a regional controllability result. We set

q(t, a, x) = q(t, a, x) e - a 0 µ(r) dr . (3.11)
Then

q satisfies        ∂ q ∂t - ∂ q ∂a -L q = 0, t 0, (a, x) ∈ (0, a b ) × Ω, ∂ q ∂ν L = 0 t 0, (a, x) ∈ (0, a b ) × ∂Ω.
(3.12)

The desired conclusion of this Proposition follows as soon as we show that there exists a constant C τ > 0 such that

a 0 0 Ω q 2 (τ, a, x) dx da C τ τ 0 a 2 a 1 ω q 2 (t, a, x) dx da dt, (3.13)
for every τ > a 1 . Indeed, using (3.13) we get

a 0 0 Ω q 2 (τ, a, x) dx da   e 2 ã 0 µ(r) dr    a 0 0 Ω q 2 (τ, a, x) dx da e 2 µ L 1 [0,a 0 ] a † 0 Ω q 2 (τ, a, x) dx da C τ τ 0 a 2 a 1 ω q 2 (t, a, x) dx da dt C τ τ 0 a 2 a 1 ω q 2 (t, a, x) dx da dt,
where C τ is a generic constant depending only on τ . We have thus shown that (3.13) implies (3.9). We can thus concentrate on the remaining part of the proof in checking (3.13).

Without loss of generality, let us assume that

τ < a 2 , τ > a 2 -a 1 , a 0 ∈ (a 1 , a 2 ). (3.14)
We set b 0 = a 2τ and we split the interval (0, a 0 ) as follows (0, a 0 ) = (0, b 0 ) ∪ (b 0 , a 1 ) ∪ (a 1 , a 0 ). (3.15) As explained in the introduction, we are going to use Proposition 3.3 along the characteristics. Before doing that, let us explain why we have divided the interval (0, a 0 ) in the above way. Basically, this division depends on the point in which the trajectory γ(s) = (τs, a + s), s ∈ [0, τ ] (or equivalently the backward characteristics staring from (τ, a)) enters the observation region (a 1 , a 2 )×(0, τ ) and exits from the same region (see Fig. 1). More precisely:

• For a ∈ (0, b 0 ), the trajectory γ(s) enters the observation region for s = a 1a. As b 0 + τ < a 2 , for s = τ, γ(s) it hits the line t = 0 without leaving the observation region (blue region in Fig. 1). • For a ∈ (b 0 , a 1 ), the trajectory γ(s) enters the observation domain for s = a 1a and exits the observation region for s = a 2a < τ (red region in Fig. 1). • For a ∈ (a 1 , a 0 ) the trajectory γ(s) starts inside the observation region but it exits the region in time s = a 2a < τ (green region in Fig. 1). Let us remark that, the choices in (3.14) are made to cover all possible scenarios. Towards the end of the proof of the proposition, we shall explain how to split the interval in other cases.

In the remaining part of the proof we give upper bounds for I Ω q 2 (τ, a, x) dx da where I is successively each one of the intervals appearing in the decomposition (3.15). Upper bound on (0, b 0 ): For a.e a ∈ (0, b 0 ), we first set w(s, x) = q(s, a + τs, x) (s ∈ (0, τ ), x ∈ Ω).

Then

w satisfies            ∂w ∂s -Lw = 0, (s, x) ∈ (0, τ ) × Ω, ∂w ∂ν L = 0, (s, x) ∈ (0, τ ) × ∂Ω, w(0, x) = q(0, τ + a, x), x ∈ Ω. (3.16)
Applying Proposition 3.3, with t 0 = 0, t 1 = τ + aa 1 and T = τ , we obtain

Ω w 2 (τ, x) dx Ω w 2 (τ + a -a 1 , x) dx c 1 e c 2 τ + a -a 1 τ +a-a 1 0 ω w 2 (s, x) dxds.
In terms of q, the above inequality writes

Ω q 2 (τ, a, x) dx c 1 e c 2 τ + a -a 1 τ +a-a 1 0 ω q 2 (s, a + τ -s, x) dxds = c 1 e c 2 τ + a -a 1 τ +a a 1 ω q 2 (τ + a -s, s, x) dxds.
Integrating with respect to a over (0, b 0 ) we obtain

b 0 0 Ω q 2 (τ, a, x) dxda c 1 e c 2 τ -a 1 b 0 0 τ +a a 1 ω q 2 (τ + a -s, s, x) dxdsda = C τ a 2 a 1 b 0 s-τ ω q 2 (τ + a -s, s, x) dxdads = C τ a 2 a 1 a 2 -s 0 ω q 2 (r, s, x) dxdrds C τ τ 0 a 2 a 1 ω
q 2 (t, a, x) dxdadt. (3.17 Upper bound on (b 0 , a 1 ): For a.e. a ∈ (b 0 , a 1 ), we define w(s, x) = q(s, a + τs, x)

(s ∈ (τ + a -a 2 , τ ), x ∈ Ω). Then w satisfies            ∂w ∂s -Lw = 0, (s, x) ∈ (τ + a -a 2 , τ ) × Ω, ∂w ∂ν L = 0, (s, x) ∈ (τ + a -a 2 , τ ) × ∂Ω, w(τ + a -a 2 , x) = q(τ + a -a 2 , a 2 , x), x ∈ Ω. (3.18)
Applying Proposition 3.3 with the choice

t 0 = τ + a -a 2 , t 1 = τ + a -a 1 and T = τ, it follows that Ω w 2 (τ, x) dx Ω w 2 (τ + a -a 1 , x) dx c 1 e c 2 a 2 -a 1 τ +a-a 1 τ +a-a 2 ω w 2 (s, x) dxds.
In terms of q, the above inequality becomes

Ω q 2 (τ, a, x) dx c 1 e c 2 a 2 -a 1 τ +a-a 1 τ +a-a 2 ω q 2 (s, a + τ -s, x) dxds = C a 2 a 1 ω z 2 (τ + a -s, s, x) dxds.
Integrating with respect to a over (b 0 , a 1 ) we get

a 1 b 0 Ω q 2 (τ, a, x) dxda C a 1 b 0 a 2 a 1 ω q 2 (τ + a -s, s, x) dxdsda = C a 2 a 1 a 1 b 0 ω q 2 (τ + a -s, s, x) dxdads = C a 2 a 1 τ +a 1 -s τ +b 0 -s ω q 2 (r, s, x) dxdrds C a 2 a 1 τ 0 ω q 2 (r, s, x) dxdrds = C τ 0 a 2 a 1 ω q 2 (t, a, x) dxdadt. (3.19)
Upper bound on (a 1 , a 0 ): For a.e. a ∈ (a 1 , a 0 ), we define

w(s, x) = q(s, a + τ -s, x) (s ∈ (τ + a -a 2 , τ ), x ∈ Ω).
Then w satisfies the system (3.18). Applying Proposition 3.3 with t 0 = τ + aa 2 and t 1 = T = τ, we have that

Ω w 2 (τ, x) dx c 1 e c 2 a 2 -a τ τ +a-a 2 ω w 2 (s, x) dxds.
In terms of q, the above inequality reads as follows

Ω q 2 (τ, a, x) dx c 1 e c 2 a 2 -a τ τ +a-a 2 ω q 2 (s, a + τ -s, x) dxds = c 1 e c 2 a 2 -a a 2
a ω q 2 (τ + as, s, x) dxds.

Integrating with respect to a over (a 1 , a 0 ) we get

a 0 a 1 Ω q 2 (τ, a, x) dxda c 1 e c 2 a 2 -a 0 a 0 a 1 a 2 a ω q 2 (τ + a -s, s, x) dxdsda = C τ a 2 a 1 s a 1 ω q 2 (τ + a -s, s, x) dxdads = C τ a 2 a 1 τ τ +a 1 -s ω q 2 (r, s, x) dxdrds C τ τ 0 a 2 a 1 ω q 2 (t, a, x) dxdadt. (3.20)
Therefore, combining (3.17), (3.19) and (3.20) we get

a 0 0 Ω q 2 (τ, a, x) dxda C τ τ 0 a 2 a 1 ω q 2 (t, a, x) dx da dt, , (3.21)
Let us explain how to split the interval (0, a 0 ) in other possible cases:

• If τ < a 2 , τ < a 2 -a 1 , a 0 ∈ (a 2 -τ, a 2 ), then we use (0, a 0 ) = (0, a 1 ) ∪ (a 1 , b 0 ) ∪ (b 0 , a 0 ).
• If τ a 2 , a 0 ∈ (a 1 , a 2 ) then we use (0, a 0 ) = (0, a 1 ) ∪ (a 1 , a 0 ). This completes the proof of the proposition.

In the next proposition, we estimate q(t, 0, x). More precisely, we prove the following: Proposition 3.5. Let us assume the hypothesis of Theorem 1.1 and let τ > a 1 and η ∈ (a 1 , τ ). Then there exists a constant C > 0 such that, for every q 0 ∈ D(A * ), the solution q of the system (3.1), satisfies

τ η Ω q 2 (t, 0, x) dxdt C τ 0 a 2 a 1 ω q 2 (t, a, x) dxdadt. (3.22)
Proof. Let q be defined as in (3.11). In particular, q satisfies (3.12). Here also we are going to use Proposition 3.3 along the characteristics. Since we want to estimate q(t, 0, x) we need to consider the trajectory γ(s) = (ts, s), s t τ (or equivalently the backward characteristic stating from (t, 0)). If τ < a 1 , the trajectory γ(s) never reaches the observation region (0, τ ) × (a 1 , a 2 ) (see Fig. 2). This is why we choose τ > a 1 . Without loss of generality, let us assume that τ > a b , η < a b and a 2 a b .

Case 1: For a.e. t ∈ (a b , τ ), we define

w(s, x) = q(s, t -s, x), s ∈ (t -a b , t), x ∈ Ω. (3.23) Then w satisfies            ∂w ∂s -Lw = 0 ((s, x) ∈ (t -a b , t) × Ω), ∂w ∂ν L = 0 ((s, x) ∈ (t -a b , t) × ∂Ω), w(t -a b , x) = q(t -a b , a b , x) (x ∈ Ω).
(3.24)

Using Proposition 3.3, with t 0 = ta b , t 1 = ta 1 and T = t, we obtain

Ω w 2 (t, x) dx Ω w 2 (t -a 1 , x) dx c 1 e c 2 a b -a 1 t-a 1 t-a b ω w 2 (s, x) dxds.
In terms of q the above inequality reads as

Ω q 2 (t, 0, x) dx c 1 e c 2 a b -a 1 t-a 1 t-a b ω q 2 (s, t -s, x) dxds = c 1 e c 2 a b -a 1 a b a 1 ω q 2 (t -s, s, x) dxds.
Integrating with respect to t over [a b , τ ] we obtain

τ a b Ω q 2 (t, 0, x) dxdt c 1 e c 2 a b -a 1 τ a b a b a 1 ω q 2 (t -s, s, x) dxdsdt = C a b a 1 τ a b ω q 2 (t -s, s, x) dxdtds = C a b a 1 τ -s a b -s ω q 2 (r, s, x) dxdrds C τ 0 a b a 1 ω q 2 (t, a, x) dxdadt. (3.25)
Case 2: For a.e t ∈ (η, a b ), we define

w(s, x) = q(s, t -s, x) (s ∈ (0, t), x ∈ Ω). (3.26) 0 a 1 a 2 t = a 1 t = τ η Figure 2.
An illustration of the estimate of q(t, 0, x). Here we have chosen a 2 = a b . Since τ > a 1 all the backward characteristics starting from (t, 0) enters the observation domain (the green region).

Then

w satisfies            ∂w ∂s -Lw = 0 ((s, x) ∈ (0, t) × Ω), ∂w ∂ν L = 0 ((s, x) ∈ (0, t) × ∂Ω), w(0, x) = q(0, t, x) (x ∈ Ω).
By applying Proposition 3.3, with t 0 = 0, t 1 = ta 1 and T = t, we obtain

Ω w 2 (t, x) dx Ω w 2 (t -a 1 , x) dx c 1 e c 2 t -a 1 t-a 1 0 ω w 2 (s, x) dxds.
This yields

Ω q 2 (t, 0, x) dx c 1 e c 2 t -a 1 t-a 1 0 ω q 2 (s, t -s, x) dxds = c 1 e c 2 t -a 1 t a 1 ω q 2 (t -s, s, x) dxds.
Integrating with respect to t over [η, a b ] we get

a b η Ω q 2 (t, 0, x) dxdt c 1 e c 2 η -a 1 a b η t a 1 ω q 2 (t -s, s, x) dxdsdt C η a b 0 t a 1 ω q 2 (t -s, s, x) dxdsdt = C a b a 1 a b s ω q 2 (t -s, s, x) dxdtds = C a b a 1 a b -s 0 ω q 2 (r, s, x) dxdrds C τ 0 a b a 1 ω q 2 (t, a, x) dxdadt. (3.27)
Combining, (3.25) and (3.27) we obtain

T η Ω q 2 (t, 0, x) dxdt C T 0 a 2 a 1 ω q 2 (t, a, x) dxdadt.
Note that, from the definition of q in (3.11), we have q(t, 0, x) = q(t, 0, x). Thus from the above estimate we clearly obtain (3.22).

3.1. Proof of the first main result. We are now in a position to prove Theorem 3.2, thus, consequently, our first main result in Theorem 1.1.

Proof of Theorem 3.2. Without loss of generality let us assume that a 1 < a †a 2 and a 2 a b .

Let us set δ = τ -(a 1 + a † -a 2 ). Let us choose ε < δ such that a 2 -ε > max{a 1 , a † -τ }.
Note that such a choice is always possible as τ > a 1 + a †a 2 a †a 2 (see Fig. 3).

a 2 -ε 0 a 1 a 2 a † t = a 1 + a † -a 2 t = τ t > a † -a t < a † -a τ = a 1 a † -τ Figure 3
. An illustration of the final time observability: For a ∈ (0, a 2ε) (blue region) the backward characteristics enters the observation domain. Thus we have the estimate (3.29). For a ∈ (a 2ε, a † ), the backward characteristics (green region) hits the line a = a † , gets renewed by the renewal condition β(a)q(t, 0, x) and then enters the observation domain (purple region). This is obtained in (3.32). Now

a † 0 Ω q 2 (τ, a, x) dxda = a 2 -ε 0 Ω q 2 (τ, a, x) dxda + a † a 2 -ε Ω q 2 (τ, a, x) dxda. (3.28)
By applying Proposition 3.4, we have

a 2 -ε 0 Ω q 2 (τ, a, x) dxda C τ τ 0 a 2 a 1 ω q 2 (t, a, x) dxdadt. (3.29)
Thus the theorem is proved as soon as we show that

a † a 2 -ε Ω q 2 (τ, a, x) dxda C τ τ 0 a 2 a 1 ω q 2 (t, a, x) dxdadt. (3.30)
Therefore in the sequel, we concentrate on proving (3.30). We recall that q is given by the formula (3.4). For a ∈ (a 2ε, a † ) we have τ + a > a † . So the expression of q in (3.4) yields

q(τ, a, •) = τ τ +a-a † π(a) π(a + t -s) e (t-s)A 2 β(a + t -s)V (s, •) ds, a ∈ (a 2 -ε, a † ). (3.31)
Integrating over Ω it is easy to verify that

Ω q 2 (τ, a, x) dx C τ τ τ +a-a † Ω V 2 (s, x) dxds.
Now integrating with respect to a over (a 2ε, a † ) we obtain

a † a 2 -ε Ω q 2 (τ, a, x) dxda C τ τ a 1 +(δ-ε) Ω q 2 (t, 0, x) dxdt. (3.32)
Finally using Proposition 3.5 to the above estimate we get (3.30). This completes the proof of the theorem.

Controls preserving positivity

An important issue in view of applications (namely in population dynamics) is to design controls such that the corresponding state trajectories join two different non-negative stationary states in some time τ , while preserving the positivity of the controlled trajectory for t ∈ [0, τ ]. This type of result has been proved in [START_REF] Hegoburu | Controllability with positivity constraints of the lotka-mckendrick system[END_REF] for the diffusion free Lotka-McKendrick system (in a uniform time) and in Lohéac, Trélat and Zuazua [START_REF] Lohéac | Minimal controllability time for the heat equation under unilateral state or control constraints[END_REF] for purely parabolic problems (in a time depending on an appropriate norm of the difference of the two stationary states). We prove below that the situation encountered in the latter case also applies to the problem considered in the present work. An essential ingredient in obtaining this type of result is proving the null controllability of the system by means of L ∞ controls and then "slowly" (s.t. positivity is preserved) driving, the initial state towards the desired target.

We first recall a classical estimate for the semigroup generated by a strictly elliptic operator with Neumann boundary conditions. Lemma 4.1. Let A 2 is defined in (2.2). Then the following holds

e tA 2 ϕ L ∞ (Ω) ϕ L ∞ (Ω) , for all t 0, ϕ ∈ L ∞ (Ω).
Proof. For the proof of this result we refer to Daners [START_REF] Daners | Heat kernel estimates for operators with boundary conditions[END_REF]Corollary 7.2] or Ouhabaz [START_REF] Ouhabaz | Analysis of heat equations on domains[END_REF]Corollary 4.10].

As a consequence of the above result, we show below that for every t 0, the restrictions to L ∞ ((0, a † ) × Ω) of the operator T t , where T is the semigroup constructed in Section 2, are bounded on L ∞ ((0, a † ) × Ω). In PDE terms, we have: Proposition 4.2. There exists a constant C > 0 such that the solution of (1.1) -(1.5) satisfies

p L ∞ ((0,τ )×(0,a † )×Ω)) C p 0 L ∞ ((0,a † )×Ω) + v L ∞ ((0,τ )×(0,a † )×Ω) , (4.1)
for every p 0 ∈ L ∞ ((0, a † ) × Ω) and v ∈ L ∞ ((0, τ ) × (0, a † ) × Ω).

Proof. Consider the operator A defined in (2.6) and the semigroup T on H = L 2 ((0, a † )×Ω) generated by A. Integrating along the characteristic lines, we have

T t p 0 =      π(a) π(a -t) e tA 2 p 0 (a -t, x) (t a), π(a)e tA 2 V (t -a, x) (t > a), (p 0 ∈ H). (4.2)
where π(a) = e -a 0 µ(r)dr and V (t, x) = a † 0 β(a)T t p 0 (a, x) da. Moreover, V (t, x) satisfies

V (t, x) = min{t,a † } 0 β(t -s)π(t -s)e tA 2 V (s, x) ds + a † min{t,a † } β(a) π(a) π(a -t)
e tA 2 p 0 (at, x)da.

From the above expression and using Lemma 4.1 we obtain

T t p 0 L ∞ ((0,a † )×Ω) C p 0 L ∞ ((0,a † )×Ω) , for all t ∈ [0, τ ],
where the constant C is independent of t. Finally using Duhamel's formula one can easily obtain (4.1).

We next show that under the assumptions in this section, the observability inequality in Theorem 3.2 can be strengthened to an inequality where the upper bound is the L 1 norm of the observation. More precisely, we have: Proposition 4.3. Under the assumption and with the notation in Theorem 3.2, for every τ > a 1 + a †a 2 , there exists k τ > 0 such that the solution q of (3.1) satisfies

a † 0 Ω q 2 (τ, a, x) dx 1 2 da k τ τ 0 a 2 a 1 ω |q(t, a, x)| dx da dt (q 0 ∈ H). (4.3)
The proof of the above proposition is similar to that of Theorem 3.2. The main idea is the same, i.e., to use observability inequality for parabolic equations along the characteristic lines. The difference is that now we want to observe the L 1 norm of q instead of L 2 norm of q. Thus we can not use Proposition 3.3. Rather we are going to use the below observability inequality for parabolic equations, which is a slight variation of a result from Fernandez-Cara and Zuazua [7, Proposition 3.2]. Proposition 4.4. Let τ > 0, 0 t 0 < τ and t 1 ∈ (t 0 , τ ]. Then for every w 0 ∈ L 2 (Ω), the solution w(s, x) of the Cauchy problem

           ∂w ∂s (s, x) -Lw(s, x) = 0, (s, x) ∈ (t 0 , T ) × Ω, ∂w ∂ν L = 0, (s, x) ∈ (t 0 , τ ) × ∂Ω, w(t 0 , x) = w 0 (x), x ∈ Ω, (4.4)
satisfies the estimate

Ω w 2 (τ, x) dx Ω w 2 (t 1 , x) dx c 1 e c 2 /(t 1 -t 0 ) t 1 t 0 ω |w(s, x)| dxds 2 , (4.5)
where the constants c 1 and c 2 depend on L, on Ω and on τ .

With the help of the above proposition we obtain: Proposition 4.5. Let us assume the hypothesis of Theorem 1.1 and let τ > a 1 , a 0 ∈ (0, a b ) and η ∈ (a 1 , τ ). Then there exists C τ > 0 such that, for every q 0 ∈ D(A 0 ), the solution q of (3.1), satisfies

a 0 0 Ω z 2 (τ, a, x) dx 1 2 da C τ τ 0 a 2 a 1 ω |z(t, a, x)| dx da dt. (4.6) (ii) τ η Ω q 2 (t, 0, x) dx 1 2 dt C τ τ 0 a 2 a 1 ω |q(t, a, x)| dx da dt. (4.7) (i)
Proof. The proof is similar to the one of Proposition 3.4 and Proposition 3.5. The only difference is that we have to use Proposition 4.4 instead of Proposition 3.3. As the procedure is completely similar we skip the details of the proof.

We can now prove Proposition 4.3.

Proof of Proposition 4.3. The proof of this theorem is similar to that of Theorem 3.2. Let δ and ε are defined as in the proof of Theorem 3.2. Then

a † 0 Ω q 2 (τ, a, x) dx 1 2 da = a 2 -ε 0 Ω q 2 (τ, a, x) dx 1 2 da + a † a 2 -ε Ω q 2 (τ, a, x) dx 1 2 da. (4.8)
On the other hand, using the expression of q in (3.4) we obtain

a † a 2 -ε Ω q 2 (τ, a, x) dx 1 2 da C τ τ a 1 +(δ-ε) Ω q 2 (t, 0, x) dx 1 2
dt. (4.9)

Finally, combining the above two estimates together with Proposition 4.5 we get (4.3).

In the following theorem we prove the null controllability of the system (1.1) -(1.5) by means of L ∞ controls. Besides the above ingredients, we use a classical duality argument, following closely the methodology in Micu, Roventa and Tucsnak [16, Proposition 2.5].

Theorem 4.6. With the notation and with the assumption in Theorem 1.1, for every τ > a 1 + a †a 2 and for every p 0 ∈ L ∞ ((0, a † ) × Ω) there exists a control v ∈ L ∞ ((0, τ ) × (a 1 , a 2 ) × ω) such that the solution p of (1.1)-(1.5) satisfies p(τ, a, x) = 0 for all a ∈ (0, a † ), x ∈ Ω.

(4.10)

Moreover, there exists a positive constant K τ such that for p 0 ∈ L ∞ ((0, a † ) × Ω) the control function v and the corresponding state trajectory p satisfy

p L ∞ ([0,τ ]×(0,a † )×Ω) + v L ∞ ([0,τ ]×(0,a † )×Ω) K τ p 0 L ∞ ((0,a † )×Ω) . (4.11)
Proof. We first remind the notation in Section 2, considering the pair (A, B), with A defined in (2.6) and B defined in (2.3), and denoting by T the C 0 semigroup on

H = L 2 ([0; a † ] × Ω) generated by A. Consider the subspace X of L 1 ([0, τ ] × [0, a † ] × Ω) defined by X = B * T * t q 0 | q 0 ∈ H . Given p 0 ∈ L ∞ ((0, a † ) × Ω) consider the linear functional F on X defined by F B * T * t q 0 = - a † 0 Ω p 0 (T * τ q 0 ) . (4.12)
Using (4.3), it follows F is well defined and that

|Fw| k τ |Ω| p 0 L ∞ ((0,a † )×Ω) w L 1 ((0,τ )×(0,a † )×Ω) (w ∈ X )
By the Hahn-Banach theorem, we can extend the linear functional F to a bounded linear functional

F on L 1 ([0, τ ] × [0, a † ] × Ω) such that |Fw| k τ |Ω| p 0 L ∞ ((0,a † )×Ω) w L 1 ((0,τ )×(0,a † )×Ω) , w ∈ L 1 ([0, τ ] × [0, a † ] × Ω) .
By the Riesz representation theorem there exists

v ∈ L ∞ ((0, τ ) × (a 1 , a 2 ) × ω) such that τ 0 a † 0 Ω v(t -σ)B * T * t q 0 + a † 0 Ω p 0 (T * τ q 0 ) = 0 (q 0 ∈ H).
From the above formula it follows that

τ 0 T τ -s Bv(s), q 0 H ds + T τ p 0 , q 0 H = 0 (q 0 ∈ H), (4.13)
which is equivalent to

T τ p 0 + τ 0 T τ -s Bv(s) ds = 0.
Thus the solution p of (1.1)-(1.5) corresponding to the control v constructed above satisfies (4.10). Moreover, we have that v L ∞ ([0,τ ]×(0,a † )×Ω) C p 0 L ∞ ((0,a † )×Ω) .

Therefore, using the above estimate and Proposition 4.2 we obtain (4.11).

Now we are in a position to prove our second main result.

Proof of Theorem 1.3.

Step 1. Let p s,I and p s,F be two non negative steady states of the system (1.1) -(1.5) and let v s,I and v s,F be the corresponding steady controls. We set

p s,r = 1 - r N p s,I + r N p s,F , v r,k = 1 - r N v s,I + r N v s,F (r = 0, 1, • • • , N ), (4.14)
where N ∈ N will be made precise later on. Using (1.9), we have that By Theorem 4.6, there exits u r ∈ L ∞ ((0, τ * ) × (a 1 , a 2 ) × ω) such that p r (τ * , a, x) = 0 (a ∈ (0, a † ), x ∈ Ω), and p r L ∞ ((0,τ * )×(0,a † )×Ω) K τ * p s,r-1 (a, x)p s,r (a, x) L ∞ ((0,a † )×Ω) , (4.17)

where the constant K τ * does not depend on r. In particular, if we set p r = p r + p s,r , v r = v r + v s,r ((t, a, x) ∈ (0, τ * ) × (0, a † ) × Ω), (4.18)

Comments and extensions

The main result in this section gives lower bounds for the controllability time in Theorem 1.1. We show, in particular, that the controllability time in Theorem 1.1 is sharp in the case a 1 = 0. More precisely, we have: Proposition 5.1. Under the assumptions of of Theorem 1.1, let τ < max{a 1 , a †a 2 }. Then there exists q 0 ∈ D(A * ) such that, the solution q of (3.1), satisfies

• q(t, a, x) = 0 for all (t, a, x) ∈ (0, τ ) × (a 1 , a 2 ) × Ω.

• a † 0 Ω q 2 (τ, a, x)dxda > 0.

Proof. Let (ϕ j ) j 1 be an orthonormal basis of L 2 (Ω) comprising of eigenvectors of the operator -A 2 and let (λ j ) j 1 be the corresponding eigenvalues. The solution q of (3.1) writes q(t, a, •) = ∞ j=1 q j (t, a)ϕ j ,

where         
∂q j ∂t -∂q j ∂a β(a)q j (t, 0) + (µ(a) + λ j )q j = 0, t 0, a ∈ (0, a † ), q j (t, a † ) = 0, t 0, q j (0, a) = q j 0 (a), a ∈ (0, a † )

(5.1) and q 0 (a, •) = ∞ j=1 q j 0 (a)ϕ j .

Integrating along the characteristic lines, we obtain the following expression of q j (t, a)

q j (t, a) =         
π j (a) π j (a + t) q j 0 (a + t) + t 0 π j (a) π j (a + ts) β(a + ts)V j (s) ds, t a †a, t t+a-a † π j (a) π j (a + ts) β(a + ts)V j (s) ds, t > a †a,

(5.2) where V j (t) = q j (t, 0) and π j (a) = exp λ j a + a 0 µ(r) dr .

Without loss of generality, let us assume that a 1 > 0.

Case 1: τ < a 1 Let us choose ā ∈ (τ, a 1 ) and ε > 0 such that (āε, ā + ε) ⊂ (τ, a 1 ). Let q 1 0 ∈ C ∞ c (0, a † ) such that q 1 0 = 0 for all a ∈ [0, τ ] ∪ [a 1 , a † ] and q 1 0 = 1 for all a ∈ (āε, ā + ε). Since β(a) = 0 for a ∈ [0, a b ), using the expression (5.2), we first have V 1 (t) = q 1 (t, 0) = 0 for all t ∈ [0, τ ]. Using this it is easy to see that

q 1 (t, a) =    π 1 (a) π 1 (a + t) if a ∈ (ā -t -ε, ā -t + ε), t ∈ [0, τ ] 0 if a a 1 , t ∈ [0, τ ].

Remark 1 . 4 .

 14 We denote by R = a † 0 β(a)e -a 0 µ(r) dr da the reproductive number. It is known that (see, for instance [2, Theorem 3.1])

  )

0 a 1 a 2 t = a 1 t = τ b 0 a 0 Figure 1 .

 2101 Figure 1. An illustration of the choice made in (3.14): Blue region corresponds to the interval (0, b 0 ), Red corresponds to the interval (b 0 , a 1 ), Green corresponds to the interval (a 1 , a 0 ).

p. 15 ) 2 .

 152 s,r (a, x) δ ( a.e. on [0, a * ] × Ω, r = 0, 1,• • • , N ). (4Step Without loss of generality let us assume that a 2 < a * . Let us fix τ * > a 1 + a †a 2 and consider the following control problem for r 1 ∂a-L p r + µ(a) p r = m v r ((t, a, x) ∈ (0, τ * ) × (0, a † ) × Ω), ∂ p r ∂ν L = 0 ((a, x) ∈ (0, a † ) × ∂Ω), p r (t, 0, x) = a † 0 β(a)p r (t, a, x) da ((t, x) ∈ (0, τ * ) × Ω),p r (0, a, x) = p s,r-1 (a, x)p s,r (a, x) ((a, x) ∈ (0, a † ) × Ω).

Acknowledgement: We would like to thank Nicolas Hegoburu for fruitful discussions which help us to improve the controllability time.

then p r satisfies the system (1.1)-(1.4) with p r (0, a, x) = p s,r-1 (a, x), p r (τ * , a, x) = p s,r (a, x) ((a, x) ∈ (0, a † ) × Ω, r = 1, 2, • • • , N ). (4. [START_REF] Tucsnak | Observation and control for operator semigroups[END_REF]) At this point we choose N sufficiently large to have p s,r-1p s,r L ∞ ((0,a † )×Ω) < δ K τ * ,

where δ > 0 is the constant appearing in (1.9). Using the above relation, together with (4.15) and (4.17), for a.e. (t, a, x) Step 3. We define

Then we can easily verify that (p, v) satisfies the system (1.1) -(1.4). Moreover, the conclusion of Theorem 1.3 holds with τ = N τ * . This completes the proof of the theorem. Now set q 0 (a, x) = q 1 0 (a)ϕ 1 (x). Then q(t, a, x) = 0 for all (t, a, x) ∈ (0, τ ) × (a 1 , a † ) × Ω. Moreover,

Case 2: a 1 τ < a †a 2 Let us choose ā ∈ (a 2 + τ, a † ) and ε > 0 such that (āε, ā + ε) ⊂ (a 2 + τ, a † ). Let q 1 0 ∈ C ∞ c (0, a † ) such that q 1 0 = 0 for all a ∈ [0, a 2 + τ ] and q 1 0 = 1 for all a ∈ (āε, ā + ε). As β(a) = 0, for a ∈ [0, a b ), using the expression (5.2) we obtain

Using the above identity and (5.2), we can easily obtain that q 1 (t, a) = 0 for all t ∈ (0, a b ), a ∈ (0, a 2 + τt), and

Next we can calculate q 1 for (t, a) ∈ (a b , 2a b ) × (0, a † ) with q 1 (a b , •) as initial data. Continuing this process, we obtain q 1 (t, a) = 0 for all t ∈ (0, τ ), a ∈ (0, a 2 + τt), and

Then we can proceed as Case 1 to conclude the proof of the proposition.

As a consequence of the above proposition, the following theorem follows easily:

Theorem 5.2. Under the assumption of Theorem 1.1, the system (1.1)-(1.5) is not null-controllable in time τ < max{a 1 , a †a 2 }.

Remark 5.3. The above theorem shows that the controllability time in Theorem 1.1 is sharp in the case a 1 = 0.

Before ending the paper, we describe some possible extensions to be considered in future work. First, still in the case a 1 = 0, it would be interesting to make precise the dependence of the control cost on a 2 . This could allow the extension to the diffusive case of the singular perturbation result obtained in [START_REF] Hegoburu | Controllability with positivity constraints of the lotka-mckendrick system[END_REF], which describes the behaviour of the control problem when a 2 → 0 (direct birth control) in the diffusion free case. Moreover, let us note that the methods in this work are easily adaptable to the case when the mortality and fertility rates depend on the spatial variable x, whereas their adaptation to time dependent mortality and fertility rates seems a more difficult question.

Other possible directions for future extensions of the results and methods in this work concern nonlinear problems (such as considering, for instance, mortality rates depending on the total populations), controllability issues for systems involving competing species or feedback control problems.