Commit e157ced0 authored by François Trigaux's avatar François Trigaux
Browse files

added twist

parent a515bd54
[Beam]
analysis_flag = 0
nkp = 51
niter = 1000
nstep = 1
ncond_pt = 2
nframe = 0
ncurv = 0
initFromFile = 0
L=10.0
\ No newline at end of file
3.999999999999999758e-11 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00
0.000000000000000000e+00 1.210592686002522039e-10 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00
0.000000000000000000e+00 0.000000000000000000e+00 1.210592686002522039e-10 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00
0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 7.048574590989939323e-09 0.000000000000000000e+00 0.000000000000000000e+00
0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 1.919999999999999987e-09 0.000000000000000000e+00
0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 7.679999999999999949e-09
# NOTE ON THIS MAKEFILE
# As the library is not installed in a default repo, i.e. /usr/lib or /usr/local/lib, you need to give it the path to the library at runtime. This is the rpath argument.
DIR=$(shell dirname `pwd`)
default:
cp -r ../../lib ./
gcc -g -fPIC -c -o main.o main.c -I../../include
gcc -g -o gebt main.o -lgebt -L./lib
clean:
rm -f main.o
rm -f gebt
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "beam.h"
int main(int argc, char* argv[])
{
int i;
double load;
// Check the input argument
if(argc<2)
{
printf("Please specify the uniform load value as arguments\n");
printf("Example: ./gebt 100\n");
exit(EXIT_FAILURE);
}
// Initializing the Beam
printf("Initializing the beam...\n");
Beam *bm = malloc(sizeof(Beam));
initBeam(bm);
printf("...Done!\n");
double *twist = malloc(sizeof(double)*bm->nelem);
for(i=0;i<bm->nmemb;i++)
{
twist[i] = 90.0 * M_PI/180.0 * i / ((double) (bm->nmemb-1.0));
}
beam_setTwist(bm,twist,M_PI/20.0);
free(twist);
// Creating an array with the uniform load for each member
printf("Creating the loads...\n");
load = atof(argv[1]);
beam_setEndLoad(bm,load,2); // 2 -> Fz
printf("...Done!\n");
char fname[32] = "output";
// Performing the analysis
printf("Performing the analysis...\n");
beam_analysis(bm);
beam_writeSolToFile(bm,fname);
printf("... Ok! \n");
return EXIT_SUCCESS;
}
# Launch the Euler-Bernoulli comparison for the beam
import numpy as np
import matplotlib.pyplot as plt
from readIni import IniReader
import subprocess as sp
ini = IniReader("Beam.ini");
nkp = ini.getValue("Beam","nkp",typeCast=int);
L = ini.getValue("Beam","L" ,typeCast=float);
nmemb = nkp-1
display = 1
x = np.linspace(0,1);
xcheck = np.linspace(0,1,nkp);
#a = np.diag(10**3 * np.array([1770, 1770, 1770, 8.16, 86.9, 215]));
#print(a);
#print(np.linalg.inv(a));
# Check "A VERIFICATION AND VALIDATION OF THE GEOMETRICALLY EXACT BEAM THEORY WITH LEGENDRE SPECTRAL FINITE ELEMENTS FOR WIND TURBINE BLADE ANALYSIS" for more details
# https://mountainscholar.org/bitstream/handle/11124/17010/Johnson_mines_0052N_10590.pdf?sequence=1
E = 200e9;
G = 79.3e9;
h = 0.5;
w = 0.25;
F = 4e6;
A = h*w;
Iy = h*h*h*w/12;
Iz = w*w*w*h/12;
J = h*w*w*w * 0.229
k = 5.0/6.0
stiff = np.diag(np.array([E*A, k*G*A, k*G*A, G*J, E*Iy, E*Iz]));
#stiff = 1e8 * np.diag(np.array([250, 92.449, 83.497, 1.498, 5.208, 1.302]));
comp = np.linalg.inv(stiff);
np.savetxt("Compliance.dat",comp);
sp.run(["./gebt", "-%e"%(F)]);
data_pts = np.loadtxt("output000.dat",max_rows=nkp);
data_mem = np.loadtxt("output000.dat",skiprows=nkp);
pts = data_pts;
mem = data_mem;
if(display):
#plt.plot(mem[:,0]+mem[:,3],mem[:,4]+mem[:,1]);
plt.plot(pts[:,0]+pts[:,3],pts[:,5]+pts[:,2],'.-k');
#plt.plot(mem[:,0]+mem[:,3],mem[:,5]+mem[:,2],'.r');
plt.xlabel('x');
plt.ylabel('w(x)');
plt.axis('scaled');
print(" Ansys\t GEBT \t Difference");
print("u_1 = %1.3f\t %1.3f\t (%+1.3f%%)"%(-1.134,pts[-1,3],(-1.134-pts[-1,3])/(-1.134)*100));
print("u_2 = %1.3f\t %1.3f\t (%+1.3f%%)"%(-1.714,pts[-1,4],(-1.714-pts[-1,4])/(-1.714)*100));
print("u_3 = %1.3f\t %1.3f\t (%+1.3f%%)"%(-3.584,pts[-1,5],(-3.584-pts[-1,5])/(-3.584)*100));
print("Done!");
if(display):
plt.show();
\ No newline at end of file
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "readINIFile.h"
#include "beam.h"
......@@ -32,8 +33,6 @@ void initBeam(Beam *bm)
section = "Beam";
SET_OR_DEFAULT(&(bm->analysis_flag),%d,fileBeam,section,"analysis_flag",2);
SET_OR_DEFAULT(&(bm->nkp),%d,fileBeam,section,"nkp",2);
SET_OR_DEFAULT(&(bm->nframe),%d,fileBeam,section,"nframe",0);// number of different cross-sectional frames
SET_OR_DEFAULT(&(bm->ncurv),%d,fileBeam,section,"ncurv",0);// number of sets of initial curvature/twists
SET_OR_DEFAULT(&(bm->L),%lf,fileBeam,section,"L",1.0);// number of sets of initial curvature/twists
SET_OR_DEFAULT(&(bm->initFromFile),%d,fileBeam,section,"initFromFile",0);// number of sets of initial curvature/twists
......@@ -49,6 +48,8 @@ void initBeam(Beam *bm)
bm->nmate = bm->nkp; // 1 material per point
bm->ncond_mb = bm->nmemb; // each member has one condition
bm->ndistrfun = bm->nmemb; // each member will have its own distributed load.
bm->nframe = bm->nmemb; // each member will have its own frame.
bm->ncurv = bm->nmemb;
bm->ncond_pt = 2; // condition for clamped-free beam
bm->ntimefun = 0;
......@@ -104,9 +105,9 @@ void initBeam(Beam *bm)
bm->member[i][1] = i+2; // number of point 2
bm->member[i][2] = i+1; // number of material 1
bm->member[i][3] = i+2; // number of material 2
bm->member[i][4] = 0; // frame number
bm->member[i][5] = 1; // number of division
bm->member[i][6] = 0; // initial curvature number
bm->member[i][4] = i+1; // frame number
bm->member[i][5] = 1; // number of division per member
bm->member[i][6] = i+1; // initial curvature number
}
// Set the clamped-free conditions
......@@ -182,7 +183,7 @@ void initBeam(Beam *bm)
for(i=0;i<bm->nframe;i++){
for(j=0;j<NDIM;j++){
for(k=0;k<NDIM;k++){
bm->frame[i][j][k]= (i==k) ? 1.0 : 0.0;
bm->frame[i][j][k]= (j==k) ? 1.0 : 0.0;
}
}
}
......@@ -223,9 +224,10 @@ void initBeam(Beam *bm)
{
for(i=0;i<bm->ncurv;i++)
{
// TODO - Allow for curvature
printf("Cannot set curvature\n");
exit(EXIT_FAILURE);
for(j=0;j<NDIM;j++)
{
bm->curvature[i][j] = 0.0;
}
}
}
......@@ -551,6 +553,25 @@ void initBeamFromGebtFile(Beam *bm)
}
// Assign the twist to each of the member by changing its frame
// twist_per_L [rad/arc length] is only required if there are subdivisions to the members
void beam_setTwist(Beam *bm, double *twist, double twist_per_L)
{
int i;
for(i=0;i<bm->nframe;i++){
bm->frame[i][0][0] = 1.0;
bm->frame[i][0][1] = 0.0;
bm->frame[i][0][2] = 0.0;
bm->frame[i][1][0] = 0.0;
bm->frame[i][1][1] = cos(twist[i]);
bm->frame[i][1][2] = -sin(twist[i]);
bm->frame[i][2][0] = 0.0;
bm->frame[i][2][1] = sin(twist[i]);
bm->frame[i][2][2] = cos(twist[i]);
bm->curvature[i][0] = twist_per_L;
}
}
void beam_writeSolToFile(Beam *bm,char* fileName)
{
......
......@@ -152,6 +152,8 @@ typedef struct
void initBeam(Beam *bm);
void beam_setTwist(Beam *bm, double *twist, double twist_per_L);
void beam_setLoads(Beam *bm, double *loads, int load_no);
void beam_setEndLoad(Beam *bm, double loads, int load_no);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment