MATLAB Answers

Plotting a 3D mode shape by revolving 2D datapoints (in the r-z-plane) of an axisymmetric geometry around the z-axis

조회 수: 34(최근 30일)
Christopher Knuth
Christopher Knuth 18 Apr 2021 23:16
편집: Christopher Knuth 22 Apr 2021 12:51
Dear Matlab community,
I would like to generate a 3D graph of a vibrational mode shape using 2D data of a solid of revolution (the cross-section of which lies in the r-z plane of a cylindrical coordinate system), by revolving the 2D data around the z-axis in the circumferential direction θ by 2π radians. I have a point cloud of data points of the meshed geometry where each corresponds to a node in a finite element model. From that model I also have the the three displacement components at each node, u (in radial direction r), v (in circumferential direction θ) and w(in axial direction z).
The undeformed geometry can be seen in the following graph (which is a simplification of the actual geometry I am solving for, but essentially includes similar geometric features, but us not as symmetric)
The red dots correspond to the nodes in the finite element model and the blue line corresponds to the extreme edge(s).
The displaced geometry (,,) of a mode at the plane where θ=0 looks as follows (after scaling it):
where the cyan line shows the edges of the undeformed geoemetry. I neglected the circumferential displacement here but it needs to be included in the 3D figure.
The displacements vary sinusoidally around the circumference, which needs to be included in the final 3D plot as:
Ideally, I would like to plot the 3D mode shape as shown in the figure below. It is taken directly from a finite element software, but I can't use it, since I apply further processing outside the software. But after finishing my calculations, I again end up with 3 displacement components at the same nodes and could plot the new ones. The colourmap c represents the magnitude of the overall displacement at each point as .
In the appendix I provide the data of the wheel in the data.mat file. The file is required to run the code below. It includes the nodal coordinates r and z of the above-shown undeformed 2D geometry and the three displacement components for the single mode. It would be sufficient to revolve the outer edges (blue lines in the first/second figure), instead of revolving all points. A good starting point is probably to try and plot the undeformed shape first, before going any further. Afterwards, the displacements could be added to each point and maybe it is possible to add a colourmap. I add a small piece of code below, which only includes my start. I have, however, already had problems coming up with something to plot the undeformed geoemtry, other than creating a 3D point cloud. I am very sorry, I can't provide you more. I'm not very experienced in plotting something in 3D in MATLAB and I am not sure how to tackle this problem, to be honest.
I am really hoping anyone could provide me some adivce on how to plot a 3D shape from the 2D data, as shown above. I am grateful for anything, I hope my description makes sense. This is the first time I am using the forum, I hope I was precise enough for you, to understand my problem. If not, please don't hesitate to tell me and I will do my best providing all neccesary information. Thank you very much in advance.
Best wishes,
%Variables from data.mat are:
%r: radial coordinate of each node
%z: axial coordinate of each node
%u: displacement in radial direction of each node
%v: displacement in circumferential direction of each node
%w: displacement in axial direction of each node
%Connect: Connectivity matrix (links nodes to individual finite elements)
clc;clear variables
%load data from data.mat
%angle theta around circumference in steps of pi/N
N = 30;
theta = 0:pi/N:2*pi;
%number of nodal diameters
n = 2;
%displacement magnitude or colourmap
c = sqrt(u.^2+v.^2+w.^2);
%re-scale displacements to max displacement
u = u./max(abs([u', v', w']));
v = v./max(abs([u', v', w']));
w = w./max(abs([u', v', w']));
%r,z,u,v,w on boundary or extreme edge nodes
rEdg = r(pEdge);
zEdg = z(pEdge);
uEdg = u(pEdge);
vEdg = v(pEdge);
wEdg = w(pEdge);
%rescale mode shape in plotby scale factor to look proper
scale = 0.1
%%%%%%%%%%%%%%%%%%%%% plots %%%%%%%%%%%%%%%%%%%%%%%%%%%
%% create 3D point cloud of undeformed geometry for each node/value of theta
%create arrays for each value of theta
%all nodes
R0 = r.*ones(1,length(theta));
Z0 = z.*ones(1,length(theta));
%edge nodes (comment out if all nodes shall be used)
% R0 = rEdg.*ones(1,length(theta))
% Z0 = zEdg.*ones(1,length(theta))
Th0 = theta.*ones(length(R0),1);
%transform cyl polar to cartesian coords
[X0,Y0,Z0] = pol2cart(Th0,R0,Z0);
axis equal
axis padded
%% create 3D point cloud of deformed geometry for each node/value of theta
%create arrays for each value of theta
%all nodes
% R = r+u*scale.*ones(1,length(theta)).*cos(n*theta);
% Z = z+w*scale.*ones(1,length(theta)).*cos(n*theta);
% Th = theta+v*scale.*sin(n*theta);
%edge nodes (comment out if all nodes shall used)
R = rEdg+uEdg*scale.*ones(1,length(theta)).*cos(n*theta);
Z = zEdg+wEdg*scale.*ones(1,length(theta)).*cos(n*theta);
Th = theta+vEdg*scale.*sin(n*theta);
%transform cyl polar to cartesian coords
[X,Y,Z] = pol2cart(Th,R,Z);
axis equal
axis padded
%% plot 3D as slices at angle theta increment
for ith = 1:length(theta)
u_theta = u*scale.*cos(n*theta(ith));
v_theta = v*scale.*sin(n*theta(ith));
w_theta = w*scale.*cos(n*theta(ith));
for i = 1:size(Connect,1)
Conn = (Connect(i,:));
re = r(Conn)+u_theta(Conn);
ze = z(Conn)+w_theta(Conn);
the = theta(ith)+v_theta(Conn);
[Xe,Ye,Ze] = pol2cart(the,re,ze);
patch(Xe,Ze,Ye,[0.5 0.5 0.5])
hold on
axis equal
axis([-0.6 0.6 -0.1 0.1 -0.6 0.6])


Matt J
Matt J 18 Apr 2021 23:53
I can't follow most of what you posted, but if you just want to generate a surface/solid of revolution, you can use cylinder()
  댓글 수: 2
Christopher Knuth
Christopher Knuth 22 Apr 2021 12:34
Thanks for trying to help me @Matt J. After some thoughts I was able to resolve the problem on my own. I did not use cylinder() as you suggested, I wasn't quite clear how to use in my case. But I found another way, which I describe in an answer below, with added code.
I never used the forum, can/shall I close the question by accepting my own answer? That seems slightly odd to me.
Best Chris

댓글을 달려면 로그인하십시오.

Christopher Knuth
Christopher Knuth 22 Apr 2021 11:35
편집: Christopher Knuth 22 Apr 2021 12:38
I have managed to plot the 3D figure from my 2D data by using the patch function. First, I create a 3D boundary mesh from my 2D data (only blue line) by revolving the points (red dots) in the 2D plane step-by-step by an angle increment. I preallocate an array with the corner points of each boundary element. Then each element can be plotted with the patch() function.
The raw 2D data must look like:
[in my case each point has a radial and axial coordinate in (r,z) + sinusoidally varying mechanical displacements (u,v,w) in 3D in (r,θ,z)]
The resulting figure for above-shown arbitrary mode of a disc is depicted below.
My code looks like the following, in case someone is interested in doing something similar. The data is specific to finite element vibration analysis, but the code may be adopted for similar problems for plotting a figure of a solid of revolution. It can be used for any type of cross-section, if one has the connectivity matrix that links the nodes (red dots) to the edge elements (blue lines), though I am showing a disc which might be easier implemented, as well. It might not be most efficient, but it works and does the job.
Still, if anyone knows a function which does that job directly or another way, I am happy for suggestions to improvements, although with my current code speed is not really a problem for the number elements/angles I require.
clc;clear variables;
%% numerical example data
%radial coordinate (undeformed) (all nodes in 2D plane)
r = [0.10,0.10,0.20,0.10,0.20,0.30,0.20,0.30,0.40,0.30,0.40,0.50,0.40,0.50,0.60,0.50,0.60,0.60]';
%axial coordinate (undeformed) (all nodes in 2D plane)
z = [-0.050,0,-0.050,0.050,4.16e-18,-0.050,0.050,0,-0.050,0.050,0,-0.050,0.050,-1.04e-18,-0.050,0.050,0,0.050]';
%circumferential coordinate (undeformed) for extension to 3D
Nangle = 40;
th = linspace(0,2*pi,Nangle);
%displacement in radial direction (all nodes in 2D plane)
u = [0;0;0.0140;0;0;0.00300;-0.0140;0;-0.0170;-0.00300;0;-0.0300;0.0170;0;-0.0330;0.0300;0;0.0330];
%displacement in circumferential direction (all nodes in 2D plane)
v = [0;0;-0.020;0;0;-0.019;0.020;0;-0.010;0.019;0;0.0020;0.010;0;0.0080;-0.0020;0;-0.0080];
%displacement in axial direction (all nodes in 2D plane)
w = [0;0;0.0420;0;0.0430;0.0700;0.0420;0.0730;0.0540;0.0700;0.0570;0;0.0540;0.00100;-0.0680;0;-0.0670;-0.0680];
%Connectivity matrix that links node numbers to edge elements (not the element inside domain)
Connect = [2,1;1,3;4,2;3,6;7,4;6,9;10,7;9,12;13,10;12,15;16,13;15,17;18,16;17,18];
%scaling factor for displacements in illustration
scale = 0.1;
%number of nodal diameters of circumferential variation (zero lines)
n = 1;
%maximum displacement (for scaling)
Dmax = max(abs([u',v',w']));
%max displacement magnitude (for colourmap)
Cmax = max(sqrt(u.^2+v.^2+w.^2));
%number of elements and nodes/element
NElem = size(Connect,1);
NNode = size(Connect,2);
%% preallocate matrix including node coordinates of each patch
for iElem = 1:NElem % sweep over each edge element
%connectivity for i-th element
Conn = Connect(iElem,:);
%elemental coordinates (re,ze) and displacements (ue,ve,we)
re = r(Conn);
ze = z(Conn);
ue = u(Conn);
ve = v(Conn);
we = w(Conn);
%rotate element in circle around circumferential coordinate
for ith = 1:Nangle-1 %sweep over each angle increment(except last because is at 2pi or 0)
%sinusiodally varying displacements of the nodes at angle i and angle i+1
ueth = [ue*cos(n*th(ith));flipud(ue)*cos(n*th(ith+1))];
veth = [ve*sin(n*th(ith));flipud(ve)*sin(n*th(ith+1))];
weth = [we*cos(n*th(ith));flipud(we)*cos(n*th(ith+1))];
%colour map/nodal displacement magnitude for boundary element and at each angle
ceth(:,ith) = sqrt(ueth.^2+veth.^2+weth.^2);
%nodal coordinate+displacement(re-scaled) for boundary element and at each angle
reth(:,ith) = [re;flipud(re)]+ueth/Dmax*scale;
theth(:,ith) = [th(ith)*ones(NNode,1);th(ith+1)*ones(NNode,1)]+veth/Dmax*scale;
zeth(:,ith) = [ze;flipud(ze)]+weth/Dmax*scale;
%transform cylindrical to cartesian coordinates
[Xe,Ye,Ze] = pol2cart(theth,reth,zeth);
%store solution for each boundary element and each angle in 3D array
Ctmp(:,iElem,1:size(ceth,2)) = ceth;
Xtmp(:,iElem,1:size(ceth,2)) = Xe;
Ytmp(:,iElem,1:size(ceth,2)) = Ye;
Ztmp(:,iElem,1:size(ceth,2)) = Ze;
%reshape 3D arrays to 2D (row: nodes of patch ,column: patches)
C = reshape(Ctmp,2*NNode,[]);
X = reshape(Xtmp,2*NNode,[]);
Y = reshape(Ytmp,2*NNode,[]);
Z = reshape(Ztmp,2*NNode,[]);
%% plot 3D figure
% patch(X,Y,Z,C,'FaceColor','interp','EdgeColor','none');
colormap jet
axis equal
xlim([-(max(r)+scale) max(r)+scale])
ylim([-(max(r)+scale) max(r)+scale])
zlim([-(max(z)+scale) max(z)+scale])
axis off





Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by