Generalized Solver for Beam Section Properties

はじめる前に
新機能一覧
Maple ワークシートの作成
Mapleワークシートを共有
Maple ウィンドウのカスタマイズ
Maple システムのカスタマイズ
コネクティビティ
Mathematics
数学
物理
プログラミング
グラフィックス
Student パッケージ
Science and Engineering
科学・エンジニアリング
アプリケーションと例題
リファレンス
システム
MapleSim
MapleSim ツールボックス
MapleSim ヘルプ
Tasks
Toolboxes
- BlockImporter
- Global Optimization
- Grid Computing
- Maple-NAG Connector
  - General information
  - Computational Routines
  - Getting Started
  - Example Worksheet
エラーメッセージガイド
マニュアル
数学アプリ

Using NAG Routines in Maple - Generalized Solver for Beam Section Properties

This document illustrates the use of the Maple-NAG Connector. In this example, NAG meshing routines are used to create a tool that can return the properties of a beam section with arbitrary shape, as defined by the user. Specifically, this routine will take in the profile of the section, including any voids in the section, as a set of , coordinates, and return:

•

the area,

•	the second moment of area (in general terms, and )

•	the position of the center of gravity , .

See figure 1.

Figure 1

These properties are essential for many structural and dynamic engineering applications. However, while most engineering references give you the formulas for a range of predefined shapes, it is very difficult to have a generalized solution for any profile.

This application uses the NAG d06bac (nag_mesh2d_bound) and d06abc (nag_mesh2d_delaunay) to generate a mesh of triangular elements over the input section profile. It then computes the area and center of gravity of each element and integrates the second moment of area of the section using the Parallel Axis (Steiner) Theorem. Given the sum of the elemental properties, the area and center of gravity for the whole section can be found.

Background

The second moment of area of a plane surface in relation to an arbitrary x or y axis in the plane is given by the relationships:

and

where and are the distances from their respective relative axes (see figure below) and is the area of an infinitesimally small element in the surface.

Since the surface cannot be mathematically defined, finding the closed-form integral is not possible. Therefore we use the NAG meshing routines to break down the surface, defined by the coordinates of the boundaries (for the outside edge and the internal voids), into small triangular elements and compute the first and second moments of area for the individual elements, using the Parallel Axis (or Steiner) theorem.

For each triangle, we can compute the area and centroid () from:

, and , and some simple geometry.

Given the centroid position in x,y coordinates from the relative axes it is then a matter of computing:

Moment of Areas:

Second Moment of Areas: ,

This is done for all N elements in the mesh:

Ax[tot] = sum(A[i]*x[S[i]], i = 1 .. N) , Ay[tot] = sum(A[i]*y[S[i]], i = 1 .. N) , A[tot] = sum(A[i], i = 1 .. N)

I[xx] = sum(A[i]*y[S[i]]^2, i = 1 .. N) , and I[yy] = sum(A[i]*x[S[i]]^2, i = 1 .. N)

Finally, the position of the Center of Gravity (CoG) for the surface can be calculated from:

x[CoG] = Ax[tot]/A[tot] , and y[CoG] = Ay[tot]/A[tot]

Initialization

Restart the Maple engine to remove any previous definitions and results.

restart;

Load the Maple packages needed for the example, including the NAG package which forms the basis of the Maple-NAG Connector.

>	with(LinearAlgebra):with(plots): with(NAG):

Conversion Routines

The NAG routines handle data matrices as packed arrays. The following routines allow the matrices to be handled in a more "human readable" form as matrices with column, row indexing. These are taken directly from the NAG documentation.

>	COORCH := proc(i,j) coorch[2(j-1) + i]; end proc: LINED := proc(i,j) lined[4(j-1) + i]; end proc: COOR := proc(i,j) coor[2(j-1) + i]; end proc: COORUS:= proc(i,j) coorus[2(j-1) + i]; end proc: EDGE:= proc(i,j) edge[3(j-1) + i]; end proc: CONN:= proc(i,j) conn[3(j-1) + i]; end proc:

Boundary Shape Function

The following function defines the shape of any lines that are not straight. The user can use the LINED(4,i) element as an identifier to tell the NAG boundary generation routine to use a specified shape function in . If LINED(4,i)=0, the line is assumed to be a straight line and this function is not referenced. The user can define any mathematical function to describe the shape for any LINED(4,i)>0 by editing the following Maple procedure definition.

fbnd := proc(i, x, y, comm)
  local ret_val, d1, d2, radius2, x0, xa, xb, y0:
  xa := comm['p']['XA']:
  xb := comm['p']['XB']:
  x0 := comm['p']['X0']:
  y0 := comm['p']['Y0']:
  ret_val := 0.0:
  if (i = 1) then
      # line 1,2,3, and 4: ellipse centred in (X0,Y0) with
      # XA and XB as coefficients
      d1 := (x - x0)/xa:
      d2 := (y - y0)/xb:
      ret_val := d1*d1 + d2*d2 - 1.0:
  elif (i = 2) then
      # line 24, 27, 33 and 38 are a circle centred in (X0,Y0)
      # with radius SQRT(RADIUS2)
      x0 := 20.5:
      y0 := 4.0:
      radius2 := 4.25:
      d1 := x - x0:
      d2 := y - y0:
      ret_val := d1*d1 + d2*d2 - radius2:
  elif (i = 3) then
      x0 := 17.0:
      y0 := 8.5:
      radius2 := 5.0:
      d1 := x - x0:
      d2 := y - y0:
      ret_val := d1*d1 + d2*d2 - radius2:
  elif (i = 4) then
      x0 := 17.0:
      y0 := 8.5:
      radius2 := 5.0:
      d1 := x - x0:
      d2 := y - y0:
      ret_val := d1*d1 + d2*d2 - radius2:
  elif (i = 5) then
      x0 := 19.5:
      y0 := 4.0:
      radius2 := 1.25:
      d1 := x - x0:
      d2 := y - y0:
      ret_val := d1*d1 + d2*d2 - radius2:
  end if:
  return(ret_val):
end proc:

comm := Nag_Comm():
comm['p'] := table(['XA'=(COORCH(1,2)-COORCH(1,4))/2.0,
                    'XB'=(COORCH(2,3)-COORCH(2,1))/2.0,
                    'X0'=(COORCH(1,2)+COORCH(1,4))/2.0,
                    'Y0'=(COORCH(2,3)+COORCH(2,1))/2.0 ]):

Input Data

The following defines the shape of the beam section that is to be meshed. Detailed descriptions of each parameter are given in the documentation for d06bac.

nvmax := 1000:
nedmx := 300:

coorch := Vector([9.5, -1, 33, 7.5, 9.5, 16, -14, 7.5, -4, 3, -2, 3, 2, 3, 4, 3, 2, 7, -2, 8, -4, 12, -2, 12, 2, 12, 4, 12, 7, 3, 9, 3, 13, 3, 16, 3, 9, 5, 12, 5, 7, 12, 10, 12, 18, 2, 21, 2, 17, 3, 20, 3, 16, 5, 20, 5, 15.5, 6, 16, 6, 18, 6, 21, 6, 16, 6.5, 18, 6.5, 18.5811, 10.0811, 21, 10.0811, 17, 10.7361, 20, 10.7361, 20.5, 12, 23, 12], datatype=float[8]):
lined := Vector([15, 1, 2, 1, 15, 2, 3, 1, 15, 3, 4, 1, 15, 4, 1, 1, 4, 6, 5, -1, 10, 10, 6, 0, 10, 7, 10, 0, 4, 8, 7, -3, 15, 14, 8, 0, 4, 13, 14, 0, 10, 9, 13, 0, 10, 12, 9, 0, 4, 11, 12, 0, 15, 5, 11, 0, 4, 16, 15, 0, 7, 19, 16, 0, 4, 20, 19, 0, 7, 17, 20, 0, 4, 18, 17, 0, 13, 22, 18, 0, 5, 21, 22, 0, 13, 15, 21, 0, 4, 24, 23, 0, 10, 24, 32, 2, 4, 31, 32, 0, 4, 34, 31, 0, 10, 34, 35, 3, 4, 36, 35, 0, 4, 40, 36, 0, 4, 39, 40, 0, 4, 38, 39, 0, 4, 37, 38, 0, 10, 37, 33, 4, 4, 30, 33, 0, 4, 29, 30, 0, 4, 27, 29, 0, 4, 28, 27, 0, 10, 26, 28, 5, 4, 25, 26, 0, 4, 23, 25, 0], datatype=integer[kernelopts('wordsize')/8]):
coorus := Vector([-2.6667, 3, -3.3333, 3, 3.3333, 3, 2.6667, 3], datatype=float[8]):
rate := Vector([0.95, 1.05, 0.95, 1.05, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], datatype=float[8]):
nlcomp := Vector([4, 10, 8, 18], datatype=integer[kernelopts('wordsize')/8]):
lcomp := Vector([1, 2, 3, 4, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 22, 21, 20, 19, 18, 17, 16, 15, 30, 29, 28, 27, 26, 25, 24, 23, 40, 39, 38, 37, 36, 35, 34, 33, 32, 31], datatype=integer[kernelopts('wordsize')/8]):
coor := Vector(2*nvmax, datatype=float[8]):
edge := Vector(3*nedmx, datatype=integer[kernelopts('wordsize')/8]):

>	nlines:=RowDimension(convert(coorch,Matrix))/2: nus:=RowDimension(convert(coorus,Matrix))/2: ncomp:=RowDimension(convert(nlcomp,Matrix)): itrace := -1:

Generating the Boundary Mesh

Calling the NAG nag_mesh2d_bound (d06bac) routine, we create the boundary mesh that will form the basis for the element mesh generation.

>	d06bac(coorch, lined, fbnd, coorus, nus, rate, nlcomp, lcomp, nvmax, nedmx, nvb, coor, nedge, edge, itrace, 'nlines' = nlines, 'ncomp' = ncomp, 'comm' = comm):

This returns two results arrays: coor and edge. coor contains the coordinates for the start and end of each line and edge contains reference numbers for each of the points that make up the boundary. In this case, the number of coordinates and edges are:

>	nvb,nedge;

(6.1)

Visualizing the Section

Using the index mapping functions defined above, EDGE contains index values for the coordinates found in COOR for the start and end of each line. Using the Maple listplot() command, you can build up a plot of all the individual lines and show them all using the display() command.

>	for i from 1 to nedge do e1:=EDGE(1,i):e2:=EDGE(2,i):e3:=EDGE(3,i): edgeplot[i]:=listplot([[COOR(1,e1),COOR(2,e1)],[COOR(1,e2),COOR(2,e2)]]): end do:

>	bplot:=display(seq(edgeplot[ii],ii=1..nedge),scaling=constrained,axes=NONE): bplot;

Generating a Delaunay-Voronoi Mesh

The following prepares the input data and calls the nag_mesh2d_delaunay (d06abc) routine to create a triangular mesh over the boundary defined above. Descriptions of the input parameters are given in the NAG user documentation for this routine.

>	nvint := 0: npropa := 1: itrace := 0:

>	conn := Vector(6015, datatype=integer[kernelopts('wordsize')/8]): weight := Vector([0], datatype=float[8]):

>	d06abc(nvb,nvint,nvmax,nedge,edge,nv,nelt,coor,conn,weight,npropa,itrace):

The routine returns two results arrays: coor and conn. coor contains the coordinates for all the vertices in the mesh and conn returns the indices for each triangular element. The number of vertices and elements are given by:

nv; nelt;

(8.1)

Visualizing the Mesh

Using the index mapping functions defined above, CONN contains index values for the coordinates found in COOR for the three vertices in each element. Using the Maple listplot() command, you can build up a plot of all the individual elements and show them all using the display() command.

>	for i from 1 to nelt do c1:=CONN(1,i): c2:=CONN(2,i): c3:=CONN(3,i): meshplot1[i]:=listplot([[COOR(1,c1),COOR(2,c1)],[COOR(1,c2),COOR(2,c2)],[COOR(1,c3),COOR(2,c3)],[COOR(1,c1),COOR(2,c1)]]): end do: display(seq(meshplot1[ii],ii=1..nelt),scaling=constrained,axes=NONE);

Calculation of the Section Properties

Looping through each triangular element, the area is computed and the distance of each center of gravity from the datum is found. From this, the total area and the overall center of gravity can be found. Then by using the Parallel Axis method, the second moment of area for each element can be found around the X-X and Y-Y relative axes that run through the center of gravity of the section. Summing all the individual moments will give the overall second moment of area, usually denoted as Ix and Yi.

#Intialize properties
Area:=0.0:
SumAx:=0.0:
SumAy:=0.0:
Ixx:=0.0:
Iyy:=0.0:

#Loop through all elements and calculate area and COGs
for i from 1 to nelt do

#Calculate the lengths of each side, s1, s2 and s3
c1:=CONN(1,i): c2:=CONN(2,i): c3:=CONN(3,i):
s1:=sqrt((COOR(1,c1)-COOR(1,c2))^2+(COOR(2,c1)-COOR(2,c2))^2):
s2:=sqrt((COOR(1,c2)-COOR(1,c3))^2+(COOR(2,c2)-COOR(2,c3))^2):
s3:=sqrt((COOR(1,c3)-COOR(1,c1))^2+(COOR(2,c3)-COOR(2,c1))^2):

#Calculate the area of the element
ss:=(s1+s2+s3)/2:
A[i]:=sqrt(ss*(ss-s1)*(ss-s2)*(ss-s3)):
Area:=Area+A[i]:

#Line equations for two bisecting lines, in the form y = m*x+c
#Bisectors of two sides (doesn't matter which)...
p1x:=(COOR(1,c1)+COOR(1,c2))/2;p1y:=(COOR(2,c1)+COOR(2,c2))/2;
p2x:=(COOR(1,c2)+COOR(1,c3))/2;p2y:=(COOR(2,c2)+COOR(2,c3))/2;

#Slope, m
m1:=(p1y-COOR(2,c3))/(p1x-COOR(1,c3));
m2:=(p2y-COOR(2,c1))/(p2x-COOR(1,c1));

#Interset, c
C1:=(COOR(2,c3)*p1x-p1y*COOR(1,c3))/(p1x-COOR(1,c3));
C2:=(COOR(2,c1)*p2x-p2y*COOR(1,c1))/(p2x-COOR(1,c1));

#Coordinates of the centre of gravity.
# NOTE: If the slope is infinity (vertical), do something about it!
if type(m1,numeric)=false then
     cogx[i]:=p1x:
     cogy[i]:=m2*cogx[i]+C2:
elif type(m2, numeric)=false then
     cogx[i]:=p2x:
     cogy[i]:=m1*cogx[i]+C1:
else
     cogx[i]:=(C2-C1)/(m1-m2);
     cogy[i]:=m1*cogx[i]+C1;
end if:

#Moment of area...
Ax[i]:=A[i]*cogx[i];
Ay[i]:=A[i]*cogy[i];

#Total moment of area...
SumAx:=SumAx+Ax[i]:
SumAy:=SumAy+Ay[i];

end do:

#Calculate the coordinates of the total section
COGx:=SumAx/Area:
COGy:=SumAy/Area:

#Finally, compute the second moment of area about the relative axes for each element and add them up
for i from 1 to nelt do

Iyy:=Iyy+A[i]*(cogx[i]-COGx)^2;
Ixx:=Ixx+A[i]*(cogy[i]-COGy)^2;

end do:

Results summary

First, the Center of Gravity...

>	COGx,COGy;display([bplot,pointplot([COGx,COGy],symbol=CROSS, color=RED)],axes=NORMAL);

Comparing this with the COG of the ellipse only. From the input data the major and minor diameters are..

>	De:=47;de:=17;

(11.1)

...but there are x,y offsets of...

>	Xoff:=-14;Yoff:=-1;

(11.2)

So comparing the two, we get...

>	COGx,COGy; De/2.0+Xoff,de/2.0+Yoff;

(11.3)

Comparing the area with the voids and without...

>	Area,evalf(Pi/4Dede);

(11.4)

Finally, comparing the Second Moments of Area for the input shape and just the ellipse. First, about the X-X relative axis...

>	Ixx,evalf((Pi(De/2)(de/2)^3)/4);

(11.5)

...then, about the Y-Y relative axis...

>	Iyy,evalf((Pi(de/2)(De/2)^3)/4);

(11.6)

Concluding Remarks

This is a relatively simple analysis that can be carried out with Maple and the NAG meshing routines, but would have been almost impossible without NAG. It highlights how easy it is to access NAG routines through Maple without developing C or FORTRAN programs, thus allowing the development and testing of the analysis in a matter of a few hours. Furthermore, once the structure of the output from these routines is understood, it very easy to visualize the results using plotting capabilities of Maple.