MBTR test rig with AMB bearings¶
This code example is based on a real magnetic bearing test rig (MBTR). In this MBTR a rotor is mounted on active magnetic bearings (AMB’s).
Configuration file for the MBTR with AMB¶
This is a Configuration file for a rotor system based on the magnetic bearing test rig (MBTR). The example of the Laval rotor (Simple laval rotor example) serves as a template and has been adapted and extended for the purposes of the MBTR rotor.
The first step is to define the rotor/shaft. This definition can be divided into three parts: Definition of the material parameters, the geometry of the rotor and the properties of the mesh:
Note
In this simulation, not the entire contour of the rotor is included in the rotor geometry. Instead, only the main radius of the rotor and some reinforcements (for transition stiffness purposes) are modeled in the rotor geometry. The other diameter steps of the rotor are later modeled as discs.
1cnfg.cnfg_rotor.name = 'MBTR - TestRig';
2% All units in SI
3cnfg.cnfg_rotor.material.name = 'steel';
4cnfg.cnfg_rotor.material.e_module = 211e9; % [N/m^2]
5cnfg.cnfg_rotor.material.density = 7860; % [kg/m^3]
6cnfg.cnfg_rotor.material.poisson = 0.3; % [-]
7% Rayleigh damping: D=alpha1*K + alpha2*M
8cnfg.cnfg_rotor.material.damping.rayleigh_alpha1= 1e-5;
9cnfg.cnfg_rotor.material.damping.rayleigh_alpha2= 10;
10
11%% Rotor Config
12% % Only main shaft radius, shaft steps as disc components
13% raW = 0.004;
14% % Format of the geometry definition: {[z, r_outer, r_inner], ...} ...
15% % without start- and end-node:
16% cnfg.cnfg_rotor.geo_nodes = {[0 0 0], [0 raW 0], [0.702 raW 0]};
17
18% Shaft with reinforcements at the bearing journals and at the disc
19raW = 0.004;
20riW = 0;
21raL = raW*1.2;
22riL = 0;
23raR = raW*1.2;
24riR = 0;
25% Format of the geometry definition: {[z, r_outer, r_inner], ...} without..
26% start- and end-node:
27cnfg.cnfg_rotor.geo_nodes = {[0 0 0], [0 raW riW], [0.094 raW riW], ...
28 [0.094 raL riL], [0.132 raL riL], [0.132 raW riW], ...
29 [0.350 raW riW], [0.350 raR riR], [0.376 raR riR], [0.376 raW riW],...
30 [0.604 raW riW], [0.604 raL riL], [0.642 raL riL], ...
31 [0.642 raW riW], [0.702 raW riW], [0.702 0 0]};
32
33% cnfg.cnfg_rotor.geo_nodes = {[0 0 0], [0 raW riW], [0.094 raW riW], ...
34% % Shaft
35% [0.094 raL riL], [0.132 raL riL], [0.132 raW riW], ...
36% % bearing bushing 1
37% [0.350 raW riW], [0.350 raR riR], [0.376 raR riR], [0.376 raW riW],...
38% % Rotor disc
39% [0.604 raW riW], [0.604 raL riL], [0.642 raL riL], ...
40% % bearing bushing 2
41% [0.642 raW riW], [0.702 raW riW], [0.702 0 0]}; % Shaft
42
43
44%% FEM Config
45cnfg.cnfg_rotor.mesh_opt.name = 'coarse mesh';
46% Number of refinement steps between d_min and d_max:
47cnfg.cnfg_rotor.mesh_opt.n_refinement = 10;
48cnfg.cnfg_rotor.mesh_opt.d_min= 0.002;
49cnfg.cnfg_rotor.mesh_opt.d_max = 0.04;
50% Definition of the element radius, if the geometry radius is not ...
51% constant in this section. Options: symmetric, mean, upper sum, lower sum:
52cnfg.cnfg_rotor.mesh_opt.approx = 'symmetric';
The inclusion of components (in this case bearings and discs) consists of the initialization of a component field in the cnfg-struct and a control variable called “count”. Each additional component (regardless of type) increases the count variable and is defined by a name, the position along the z-axis of the rotor, the component type and additional parameters depending on the component:
Note
The radius and the width of the disc components is only required for visualization purposes. The mass and the mass moments of inertia are declared separatly and origin from external CAD-files.
1%% Initiation of the components section in the struct
2count = 0;
3cnfg.cnfg_component = [];
4
5%% Disc
6count = count+1;
7cnfg.cnfg_component(count).name = 'AMB_rotor_L'; % incl. clamping sleeve
8cnfg.cnfg_component(count).type = 'Disc';
9cnfg.cnfg_component(count).position = 113e-3; % [m]
10cnfg.cnfg_component(count).radius = 46e-3; % [m] only for visualization
11cnfg.cnfg_component(count).width = 0e-3; % only for visualization
12cnfg.cnfg_component(count).color = AMrotorTools.TUMColors.TUMBlue; % color
13cnfg.cnfg_component(count).m = 0.85074; % disc mass [kg]
14cnfg.cnfg_component(count).Jx = 4.840066e-4; % disc mom. of ...
15 % inertia [kg*m^2]
16cnfg.cnfg_component(count).Jz = 4.036839e-4; % disc mom. of ...
17 % inertia [kg*m^2]
18cnfg.cnfg_component(count).Jp = cnfg.cnfg_component(count).Jz; % disc ...
19 %polar mom. of inertia [kg*m^2]
20
21count = count+1;
22cnfg.cnfg_component(count).name = 'Center disc'; % incl. clamping ...
23 % sleeve, screws, nuts,..
24cnfg.cnfg_component(count).type = 'Disc';
25cnfg.cnfg_component(count).position = 363e-3; % [m]
26cnfg.cnfg_component(count).radius = 96e-3; % [m], only for visualization
27cnfg.cnfg_component(count).width = 10e-3; % only for visualization
28cnfg.cnfg_component(count).color = AMrotorTools.TUMColors.TUMBlue; % color
29cnfg.cnfg_component(count).m = 1.276499; % disc mass [kg]
30cnfg.cnfg_component(count).Jx = 1.487212e-3; % disc mom. of ...
31 % inertia [kg*m^2]
32cnfg.cnfg_component(count).Jz = 2.832751e-3; % disc mom. of ...
33 % inertia [kg*m^2]
34cnfg.cnfg_component(count).Jp = cnfg.cnfg_component(count).Jz; % disc ...
35 %polar mom. of inertia [kg*m^2]
36
37count = count+1;
38cnfg.cnfg_component(count).name = 'AMB_rotor_R'; % incl. clamping sleeve
39cnfg.cnfg_component(count).type = 'Disc';
40cnfg.cnfg_component(count).position = 623e-3; % [m]
41cnfg.cnfg_component(count).radius = 46e-3; % [m], only for visualization
42cnfg.cnfg_component(count).width = 0e-3; % only for visualization
43cnfg.cnfg_component(count).color = AMrotorTools.TUMColors.TUMBlue; % color
44cnfg.cnfg_component(count).m = 0.85074; % disc mass [kg]
45cnfg.cnfg_component(count).Jx = 4.840066e-4; % disc mom. of ...
46 % inertia [kg*m^2]
47cnfg.cnfg_component(count).Jz = 4.036839e-4; % disc mom. of ...
48 % inertia [kg*m^2]
49cnfg.cnfg_component(count).Jp = cnfg.cnfg_component(count).Jz; % disc ...
50 %polar mom. of inertia [kg*m^2]
51
52%% Bearings
53% for in MBTR integrated rotor
54
55kSchaetzung = 3.5e4; % estimated stiffness [N/m]
56dSchaetzung = 90; % estimated damping [Ns/m]
57kML = kSchaetzung;
58dML = dSchaetzung;
59
60count = count +1;
61cnfg.cnfg_component(count).name = 'AxBearingLeft';
62cnfg.cnfg_component(count).position=113e-3; % [m]
63cnfg.cnfg_component(count).type='Bearings';
64cnfg.cnfg_component(count).subtype='SimpleAxialBearing';
65cnfg.cnfg_component(count).stiffness=1e10; % [N/m]
66cnfg.cnfg_component(count).damping = 0; % [Ns/m]
67
68count = count + 1;
69cnfg.cnfg_component(count).name = 'TorqueBearingLeft';
70cnfg.cnfg_component(count).position=113e-3; % [m]
71cnfg.cnfg_component(count).type='Bearings';
72cnfg.cnfg_component(count).subtype='SimpleTorqueBearing';
73cnfg.cnfg_component(count).stiffness=1e10; % [N/m]
74cnfg.cnfg_component(count).damping = 0; % [Ns/m]
The inclusion of sensors consists of the initialization of a sensor field in the cnfg-struct and a control variable called “count”. Each additional sensor (regardless of type) increases the count variable and is defined by a name, the position along the z-axis of the rotor, the sensor type and additional parameters depending on the sensor:
1%% Initialization of the sensor section in the struct
2cnfg.cnfg_sensor=[];
3count = 0;
4
5count = count + 1;
6cnfg.cnfg_sensor(count).name = 'EddyL1';
7cnfg.cnfg_sensor(count).position=0.071; % [m]
8cnfg.cnfg_sensor(count).type='Displacementsensor';
9
10count = count + 1;
11cnfg.cnfg_sensor(count).name = 'DisplAMBL';
12cnfg.cnfg_sensor(count).position=0.113; % [m]
13cnfg.cnfg_sensor(count).type='Displacementsensor';
14
15count = count + 1;
16cnfg.cnfg_sensor(count).name = 'EddyL2';
17cnfg.cnfg_sensor(count).position=0.155; % [m]
18cnfg.cnfg_sensor(count).type='Displacementsensor';
19
20count = count + 1;
21cnfg.cnfg_sensor(count).name='Laser';
22cnfg.cnfg_sensor(count).position=0.363; % [m]
23cnfg.cnfg_sensor(count).type='Displacementsensor';
24
25count = count + 1;
26cnfg.cnfg_sensor(count).name = 'EddyR1';
27cnfg.cnfg_sensor(count).position=0.581; % [m]
28cnfg.cnfg_sensor(count).type='Displacementsensor';
29
30count = count + 1;
31cnfg.cnfg_sensor(count).name='DisplAMBL';
32cnfg.cnfg_sensor(count).position=0.623; % [m]
33cnfg.cnfg_sensor(count).type='Displacementsensor';
34
35count = count + 1;
36cnfg.cnfg_sensor(count).name = 'EddyR2';
37cnfg.cnfg_sensor(count).position=0.665; % [m]
38cnfg.cnfg_sensor(count).type='Displacementsensor';
39
40count = count + 1;
41cnfg.cnfg_sensor(count).name='ControllerForceAMBL';
42cnfg.cnfg_sensor(count).position=0.113; % [m]
43cnfg.cnfg_sensor(count).type='ControllerForceSensor';
44
45count = count + 1;
46cnfg.cnfg_sensor(count).name='ControllerForceAMBR';
47cnfg.cnfg_sensor(count).position=0.623; % [m]
48cnfg.cnfg_sensor(count).type='ControllerForceSensor';
49
50count = count + 1;
51cnfg.cnfg_sensor(count).name='ExcitationForceAMBL';
52cnfg.cnfg_sensor(count).position=0.113; % [m]
53cnfg.cnfg_sensor(count).type='ForceLoadPostSensor';
54
55count = count + 1;
56cnfg.cnfg_sensor(count).name='ExcitationForceAMBR';
57cnfg.cnfg_sensor(count).position=0.623; % [m]
58cnfg.cnfg_sensor(count).type='ForceLoadPostSensor';
The inclusion of loads consists of the initialization of a load field in the cnfg-struct and a control variable called “count”. Each additional load (regardless of type) increases the count variable and is defined by a name, the position along the z-axis of the rotor, the load type and additional parameters depending on the load:
1%% Initialization of the load section in the struct
2cnfg.cnfg_load=[];
3count = 0;
4
5% Look up table force
6count = count + 1;
7cnfg.cnfg_load(count).name='Random force AMBL';
8cnfg.cnfg_load(count).position=0.113; % [m]
9cnfg.cnfg_load(count).LUT.time = 0:1e-3:5; % time range [s]
10cnfg.cnfg_load(count).LUT.Fx = rand(length(cnfg.cnfg_load(count) ...
11 .LUT.time),1); % trans. force x [N]
12cnfg.cnfg_load(count).LUT.Fy = zeros(length(cnfg.cnfg_load(count) ...
13 .LUT.time),1); % trans. force y [N]
14cnfg.cnfg_load(count).LUT.Fz = zeros(length(cnfg.cnfg_load(count) ...
15 .LUT.time),1); % trans. force z [N]
16cnfg.cnfg_load(count).LUT.Mx = zeros(length(cnfg.cnfg_load(count) ...
17 .LUT.time),1); % torque x [Nm]
18cnfg.cnfg_load(count).LUT.My = zeros(length(cnfg.cnfg_load(count) ...
19 .LUT.time),1); % torque y [Nm]
20cnfg.cnfg_load(count).LUT.Mz = zeros(length(cnfg.cnfg_load(count) ...
21 .LUT.time),1); % torque z [Nm]
22cnfg.cnfg_load(count).type='Force_timevariant_LUT';
23
24% Look up table force
25count = count + 1;
26cnfg.cnfg_load(count).name='Random force AMBR';
27cnfg.cnfg_load(count).position=0.623; % [m]
28cnfg.cnfg_load(count).LUT.time = 0:1e-3:5;
29cnfg.cnfg_load(count).LUT.Fx = rand(length(cnfg.cnfg_load(count) ...
30 .LUT.time),1); % trans. force x [N]
31cnfg.cnfg_load(count).LUT.Fy = zeros(length(cnfg.cnfg_load(count) ...
32 .LUT.time),1); % trans. force y [N]
33cnfg.cnfg_load(count).LUT.Fz = zeros(length(cnfg.cnfg_load(count) ...
34 .LUT.time),1); % trans. force z [N]
35cnfg.cnfg_load(count).LUT.Mx = zeros(length(cnfg.cnfg_load(count) ...
36 .LUT.time),1); % torque x [Nm]
37cnfg.cnfg_load(count).LUT.My = zeros(length(cnfg.cnfg_load(count) ...
38 .LUT.time),1); % torque y [Nm]
39cnfg.cnfg_load(count).LUT.Mz = zeros(length(cnfg.cnfg_load(count) ...
40 .LUT.time),1); % torque z [Nm]
41cnfg.cnfg_load(count).type='Force_timevariant_LUT';
The pidController is defined together with the definition of the AMB (see the following codeblock about the AMB’s). Nevertheless the initialization of the pidController is necessary to avoid errors. The initialization consists of the assignment of a pidController field in the cnfg-struct and the definition of a control variable “count”:
1%% Initialization of the pid-controller section in the struct
2cnfg.cnfg_pid_controller=[];
3count = 0;
The inclusion of active magnetic bearing (AMB) consists of the initialization of a active magnetic bearing (AMB) field in the cnfg-struct and a control variable called “count”. Each additional active magnetic bearing (AMB) (regardless of type) increases the count variable and is defined by a name, the position along the z-axis of the rotor and additional parameters from the active magnetic bearing (AMB):
1%% Initialization of the active magnetic bearing section in the struct
2cnfg.cnfg_activeMagneticBearing = [];
3count = 0;
4
5electricalP = 5000; % proportional controller parameter [A/m]
6electricalI = 1500; % integral controller parameter [A/ms]
7electricalD = 5; % derivative controller parameter [As/m]
8
9count = count + 1;
10cnfg.cnfg_activeMagneticBearing(count).name = 'AMB1';
11cnfg.cnfg_activeMagneticBearing(count).position = 0.113; % [m]
12cnfg.cnfg_activeMagneticBearing(count).pidType = 'pidControllerLinear'; % .
13 % .. controller type
14cnfg.cnfg_activeMagneticBearing(count).kx = -1e5; % stiffness factor ...
15 % displacement [N/m]
16cnfg.cnfg_activeMagneticBearing(count).ki = 50; % stiffness factor ...
17 % displacement [A/N]
18cnfg.cnfg_activeMagneticBearing(count).targetDisplacementX = 0; % [m]
19cnfg.cnfg_activeMagneticBearing(count).targetDisplacementY = 0; % [m]
20cnfg.cnfg_activeMagneticBearing(count).electricalP = electricalP; % [A/m]
21cnfg.cnfg_activeMagneticBearing(count).electricalI = electricalI; % [A/ms]
22cnfg.cnfg_activeMagneticBearing(count).electricalD = electricalD; % [As/m]
23
24count = count + 1;
25cnfg.cnfg_activeMagneticBearing(count).name = 'AMB2';
26cnfg.cnfg_activeMagneticBearing(count).position = 0.623;
27cnfg.cnfg_activeMagneticBearing(count).pidType = 'pidControllerLinear'; % .
28 % .. controller type
29cnfg.cnfg_activeMagneticBearing(count).kx = -1e5; % stiffness factor ...
30 % displacement [N/m]
31cnfg.cnfg_activeMagneticBearing(count).ki = 50; % stiffness factor ...
32 % displacement [A/N]
33cnfg.cnfg_activeMagneticBearing(count).targetDisplacementX = 0; % [m]
34cnfg.cnfg_activeMagneticBearing(count).targetDisplacementY = 0; % [m]
35cnfg.cnfg_activeMagneticBearing(count).electricalP = electricalP; % [A/m]
36cnfg.cnfg_activeMagneticBearing(count).electricalI = electricalI; % [A/ms]
37cnfg.cnfg_activeMagneticBearing(count).electricalD = electricalD; % [As/m]
Simulation file for the MBTR with AMB¶
This is a Simulation file for an investigation of a real AMB test rig (MBTR). Of interest is the Run-up time behaviour (from 0 to 1000 rpm) under the influence of external random excitations at the AMB’s.
Closing of all previous figures and cleaning of the workspace:
1close all
2clear all
3% clc
Import of the file path and the corresponding Configuration file:
1%% Import and formating of the figures
2
3import AMrotorSIM.* % path
4Config_Sim_MLPS % corresponding cnfg-file
5Janitor = AMrotorTools.PlotJanitor(); % Instantiation of class PlotJanitor
6Janitor.setLayout(2,3); %Setting layout of the figures
Assembly and visualization of the model:
Note
The assembly of the model is the most important part of the Simulation file and must be done before the analyses.
1%% Assembly of the rotordynamic model
2
3r=Rotorsystem(cnfg, ...
4'MBTR-Rotor with PID-Controller and negative stiffness for the AMBs');
5 % Instantiation of class Rotorsystem
6r.assemble; % Assembly of the model parts, considering the ...
7 % components (sensors,..) from the cnfg-file
8r.rotor.assemble_fem; % Assembly of the global system matrices: M, D, G, K
9
10%% Visualization of the assembled rotor model
11
12r.show; % lists the associated components of the model in teh Matlab ...
13 % Command Window
14
15r.rotor.show_2D(); % Plot of a side view of the rotor elements
16
17g=Graphs.Visu_Rotorsystem(r); % Instantiation of class Visu_Rotorsystem
18g.show(); % Plot of a 3D-isometry of the rotor with sensors, loads,...
2D side view of the rotor (left) and 3D isometry (right):
Note
The diameter steps in the rotor geometry (in the 2D side view) represent the reinforcements (check Configuration file for the MBTR) of the main rotor and not the disc component as seen in the 3D isometry.
Execution of the Run-up time integration (from 0 to 1000 rpm for 0.2 seconds) with external excitation (random on both AMB’s only in x-direction) with visualization:
1%% Running Time Simulation
2%% Stationary and runup with avaliable calculation methods and visualization
3% St_Lsg = Experiments.Stationaere_Lsg(r,0,(0:0.001:1)); % In...
4 %stantiation of class Stationaere_Lsg (Stationary solution)
5
6% St_Lsg.compute_ode15s_ss; % ode15s - method
7% St_Lsg.compute_newmark; % newmark - method
8
9Runup = Experiments.Hochlaufanalyse(r,[0,1e3],(0:0.001:0.2)); % In...
10 %stantiation of class Hochlaufanalyse (Runup)
11Runup.compute_ode15s_ss % ode15s - method
12
13%% Export and visualization of the results
14%% Export
15
16d = Dataoutput.TimeDataOutput(Runup); % Instantiation of class ...
17 % TimeDataOutput
18
19% dataset = d.compose_data(); % container: rpm ->
20 % (n,t,allsensorsxy)
21dataset = d.compose_data_sensor_wise();
22datasetName = 'MLPS_pid_5s';
23d.save_data(dataset,datasetName);
24struct = d.convert_data_to_struct_sensor_wise(dataset);
25d.save_data(struct,datasetName);
26d.write_data_to_unv(datasetName);
27
28%% Visualizing results
29
30%Lsg = St_Lsg;
31Lsg = Runup;
32
33t = Graphs.TimeSignal(r, Lsg); % Instantiation of class TimeSignal
34f = Graphs.Fourierdarstellung(r, Lsg); % Instantiation of class ...
35 % Fourierdarstellung
36 for sensor = r.sensors
37 t.plot(sensor); % Time signal
38 f.plot(sensor); % Fourier
39 Janitor.cleanFigures(); % Formatting of the figures
40 end
41
Time signal (x-direction) of the eddy-current sensor on the left AMB (left) and the corresponding Fourier analysis (right):
