von Mises Stress due to Centrifugal force on a rotor

Hallo,

I am trying to calculate and visualise von Mises Stress due to the centrifugal forces on a rotor using Finite Element Method.

I have modelled the rotor geometry and specified the fixed boundary condition. How do I specify the centrifugal force using structuralBoundaryLoad fucntion? Is there any other convenient way?

Thanks in advance!

Sridhar

Answer 1

[Ravi Kumar]

Hi Sridhar,

Centrifugal load is not yet supported in structural workflow. However, for a specific geometry you can define the centrifugal load as f-coefficient in the generic PDE workflow. Following example solves a spinning disk using generic workflow. Updated the function spinEffect as per your need.

% Create a model, assign a disk geometry, and generate mesh.
model = createpde(3);
gm = multicylinder([0.01,0.05],1E-3,'Void',[1,0]);
model.Geometry = gm;
generateMesh(model);

% Plot the geometry and note the face number of the hub to apply fixed BC.
pdegplot(model,'FaceLabels','on')
applyBoundaryCondition(model,'dirichlet','Face',3,'u',[0,0,0]);

% Specify material properties that defines the PDE as coefficients.
E = 210E9;
nu = 0.3;
rho = 7800;
c = elasticityC3D(E,nu);

% Specify the c-coefficient and centrifugal body load as f-coefficient.
specifyCoefficients(model,'m',0,'d',0,'c',c,'a',0,'f',@spinEffect);

% Solve the model and plot resultant displacement
R = solvepde(model);
figure
title('Displacement of spinning disk')
pdeplot3D(model,'ColorMapData',sqrt(R.NodalSolution(:,1).^2 + R.NodalSolution(:,2).^2+R.NodalSolution(:,3).^2))

% Find stress components using evaluateCGradient function. 
[cgradx,cgrady,cgradz] = evaluateCGradient(R);

% Note that each output argument from evaluateCGradient is a matrix of size
% number of nodes x 3, sxx = cgradx(:,1), sxy = cgradx(:,2)... and so
% on.

sxx = cgradx(:,1);
sxy = cgradx(:,2);
sxz = cgradx(:,3);

syx = cgrady(:,1);
syy = cgrady(:,2);
syz = cgrady(:,3);

szx = cgradz(:,1);
szy = cgradz(:,2);
szz = cgradz(:,3);

sVonMises = sqrt( 0.5*( (sxx-syy).^2 + (syy -szz).^2 +...
           (szz-sxx).^2) + 3*(sxy.^2 + syz.^2 + szx.^2));

figure
title('von Mises stress of spinning disk')
pdeplot3D(model,'ColorMapData',sVonMises)

function f = spinEffect(region,~)
xcg = 0; ycg = 0;
[theta,r] = cart2pol(region.x-xcg,region.y-ycg);
rho = 7800; % density.
omega = 1047; % rad/sec, about 10,000 rpm.
f = rho*omega^2.*[r.*cos(theta);r.*sin(theta);zeros(1,numel(region.y))];
end

Regards, Ravi

[Sridhar Balasubramanian]

Thanks a lot Ravi!

really helpful!

Sridhar.

[Bo Zhang]

Hi Ravi,

Is it possible to perform a similar analysis in 2D?

Thanks,

Bo

[Ravi Kumar]

Certainly. Refer to one of the 2-D examples on how to create the required geometry and change the c-coefficient for plane-stress: c = [2*G+mu;0;G;0;G;mu;0;G;0;2*G+mu];

You can easily modify the above example by eliminating z-dependency.

[Bo Zhang]

Hi Ravi, Thank you for the comment. Could you explain how you get the c coefficient? I went through the matlab c Coefficient page, but I think the c vector should have 16 elements for 2D 2 PDE equations (4*2^2). Why does your vector only contains 10 elements?

[Ravi Kumar]

Hi Bo,

For an isotropic material, the symmetric tensor can be expressed compactly by specifying only symmetric portion, which corresponds to 2N(2N+1)/2 = 10 for N = 2, format in the c-Coefficient's documentation page.

Regards,

Ravi

[Ahmed MALIM]

Hi Ravi,

I want to adapt your code ( for the centrifugal force) for my case which is a rotating Propeller

would you please explain the @spinEffect Function:

what deos mean xcg = 0; ycg = 0; and [theta,r] = cart2pol(region.x-xcg,region.y-ycg);

and how can I modify f expression if the axis of rotation is X or Y ( my propeller rotates around X axis and it is shifted from the origine (0, 0, 0) to X =~ 0.265 )

You find in attachment my propeller fig

Best regards,

Ahmed

[Ravi Kumar]

If your model's center of gravity is not at (0,0) you can use xcg and ycg to offset so you can use the cart2pool. In your case xcg = 0.265 and ycg =0, I think.

Regards,

Ravi

[Ahmed MALIM]

Thanks

What about the other axis of rotation and why did you put the z component of the centrifugal force equal to zero ??

f = rho*omega^2.*[r.*cos(theta);r.*sin(theta);zeros(1,numel(region.y))];

the expression [r.*cos(theta);r.*sin(theta);zeros(1,numel(region.y))] should present the distance between each element or node and the axis of rotation? isn't it ? would you please explain it to me?

Best regards,

Ahmed

[Lukas]

Hello, I want to simulate the same, but I want to have magnets in it additionally. How would you solve this? Two different materials in one geometry. Would you separate the rotor model in areas of magnet and iron, and just solve two different problems. Or is there a way to calculate it as a whole rotor.

Best regards,

Lukas