- Through Vectorize Operations: Where possible, use vectorized operations instead of loops. This often speeds up execution in MATLAB.
- By Combining “StructuralIC” Calls: We can Minimize the number of calls to “structuralIC” by grouping the nodes and applying velocities in a single call if feasible.
When defining StructuralIC (initial conditions in a transient-solid), is there a way to give each node a different velocity without solving takinhg too much time?
3 views (last 30 days)
Show older comments
I have the following code where to one of the gears in the two-spur-gear model has angular velocity 10rpm. To define angular velocity, I have given each point on the gear a unique instantaneous velocity as such:
angular_velocity=(10/60)*2*pi;
model=createpde('structural','transient-solid');
importGeometry(model,"full_50_teeth_gear_assembly2_3d_slight_offset.stl");
[m,n]=size(model.Geometry.Vertices);
% recognizing which geometric points belong to gear 1
t=0;
B=zeros(m/2,n);
index_B_used=zeros(m/2,1);
for i=1:m
if model.Geometry.Vertices(i,3)==11 || model.Geometry.Vertices(i,3)==-9
t=t+1;
B(t,:)=model.Geometry.Vertices(i,:);
index_B_used(t,1)=i;
else
end
end
% find centre of gear 1
g1_middle_x=(min(B(:,1))+max(B(:,1)))/2;
g1_middle_y=(min(B(:,2))+max(B(:,2)))/2;
[lb1,lb2]=size(B);
%define gear 2 to have 0 initial velocity and whole model to have 0 initial displacement
structuralIC(model,"Cell",2,"Velocity",[0;0;0]);
structuralIC(model,"Displacement",[0;0;0]);
%give each gear 1 point a unique instantaneous velocity
v_centre_to_point=zeros(3,lb1);
for k=1:lb1
v_centre_to_point(1:2,k)=[(B(k,1)-g1_middle_x);(B(k,2)-g1_middle_y)];
v_centre_to_point(3,k)=0;
end
structuralIC(model,"Vertex",index_B_used(1:lb1,1),"Velocity",[0;0;0]);
for h=1:lb1
structuralIC(model,"Vertex",index_B_used(h,1),"Velocity",angular_velocity.*v_centre_to_point(:,h));
end
%then mesh and constrain
mesh=generateMesh(model,'Hmax',100,GeometricOrder="linear");
%material properties
structuralProperties(model,'Cell',1:2,'YoungsModulus',200e9,'PoissonsRatio',0.29,'MassDensity',7870);
structuralBC(model,'Face',1,'Constraint','roller');
structuralBC(model,'Face',204,'Constraint','roller');
%add resistive toque to gear 2
gear_2_torque=200;
structuralBoundaryLoad(model,'Face',204,'SurfaceTraction',[0;0;-gear_2_torque]);
%solve
Q = meshQuality(mesh);
A=min(Q);
tic
tlist = linspace(0,1,2);
results = solve(model,tlist);
figure(6)
pdeplot3D(model,'ColorMapData',results.VonMisesStress)
title('von Mises stress')
toc
fprintf('Minimum quality of mesh (1=perfect) = %f \n',A)
0 Comments
Answers (1)
Karan Singh
on 18 Jun 2024
Hi Tanishq,
From what I can think off and according to my knowledge, we can optimize the code in 2 ways:
The code would be modified as follows:-
% Define constants
angular_velocity = (10/60) * 2 * pi;
% Create the structural model
model = createpde('structural', 'transient-solid');
% Import geometry
importGeometry(model, "full_50_teeth_gear_assembly2_3d_slight_offset.stl");
% Get geometry vertices
vertices = model.Geometry.Vertices;
[m, n] = size(vertices);
% Identify gear 1 vertices
gear1_mask = (vertices(:, 3) == 11 | vertices(:, 3) == -9);
gear1_vertices = vertices(gear1_mask, :);
index_B_used = find(gear1_mask);
% Calculate the center of gear 1
g1_middle_x = mean(gear1_vertices(:, 1));
g1_middle_y = mean(gear1_vertices(:, 2));
% Define initial conditions for gear 2 and the entire model
structuralIC(model, "Cell", 2, "Velocity", [0; 0; 0]);
structuralIC(model, "Displacement", [0; 0; 0]);
% Calculate unique instantaneous velocities for gear 1 vertices
v_centre_to_point = [gear1_vertices(:, 1) - g1_middle_x, gear1_vertices(:, 2) - g1_middle_y, zeros(size(gear1_vertices, 1), 1)]';
initial_velocities = angular_velocity .* v_centre_to_point;
% Apply initial velocities to gear 1 vertices
structuralIC(model, "Vertex", index_B_used, "Velocity", initial_velocities);
% Generate the mesh
mesh = generateMesh(model, 'Hmax', 100, GeometricOrder="linear");
% Define material properties
structuralProperties(model, 'Cell', 1:2, 'YoungsModulus', 200e9, 'PoissonsRatio', 0.29, 'MassDensity', 7870);
% Apply boundary conditions
structuralBC(model, 'Face', 1, 'Constraint', 'roller');
structuralBC(model, 'Face', 204, 'Constraint', 'roller');
% Apply a resistive torque to gear 2
gear_2_torque = 200;
structuralBoundaryLoad(model, 'Face', 204, 'SurfaceTraction', [0; 0; -gear_2_torque]);
% Solve the model
Q = meshQuality(mesh);
A = min(Q);
tic;
tlist = linspace(0, 1, 2);
results = solve(model, tlist);
% Plot the results
figure(6);
pdeplot3D(model, 'ColorMapData', results.VonMisesStress);
title('von Mises stress');
toc;
fprintf('Minimum quality of mesh (1=perfect) = %f \n', A);
Here I have Calculated “v_centre_to_point” and “initial_velocities” using vectorized operations instead of a loop and applied all initial velocities for gear 1 in a single call to “structuralIC” by leveraging the vectorized “initial_velocities”.
As I don’t have the required STL file, I can’t verify but vectorization is better than loop, while coding in MATLAB. You can read more about the same here - https://www.mathworks.com/matlabcentral/answers/2032009-what-is-vectorization-why-use-vectorization-instead-of-loops
Hope it helps!
0 Comments
See Also
Categories
Find more on Geometry and Mesh in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!