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?

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)

Answer 1

[Karan Singh]

Hi Tanishq,

From what I can think off and according to my knowledge, we can optimize the code in 2 ways:

1. Through Vectorize Operations: Where possible, use vectorized operations instead of loops. This often speeds up execution in MATLAB.
2. 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.

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!