Good morning
I’m trying to build a simple model of a disc, constrained at the center, with a dynamic axial load applied at its radius. Ultimately, I’d like to use the displacement at the outer diameter as a chekc on resonant frequency.
However, I am stuck with an error I do not understand with the solve function.
Any advice would be greatly appreciated.
THANK YOU!
here is the code:
/* Forced vibration of a disc to look at dynamic transmissibility.
Axisymmetric modeling to try to identify resonance to match with predictions
*/
/*
------ Frequency Definitions ------
*/
real[int] ff=[10.0,280.,550.,650.,800.]; // the table of selected frequencies
// include “freq_disc.edp”;
/*
------ Geometry -------
*/
real a=4.752.54/100; // disc radius
real l=0.12.54/100; // disc thickness
/*
------ Output Filename -------
*/
string ftitle = “Point_transmissibility.dat”;
/*
------ Material Definitions and Loading-------
*/
complex V0=0.0, V1=1.0; // electrode potentials
real rho = 7850; // material density
complex E = 210E9; // Young’s modulus steel
complex nu = 1/3; // Poisson ratio steel
// Stiff - Stiffness Matrix
func Stiff = (E/((1+nu)(1-2nu)))* [[1-nu, nu, nu, 0 ],
[nu, 1-nu, nu, 0 ],
[nu, nu, 1-nu, 0 ],
[0 , 0 , 0 , 1-2*nu]];
// force applied at the end of the disc
complex force = 1.0;
/*
------ Mesh Definition and Creation -------
*/
int MM=10, NN=5; // grid resolution
mesh Sp = square(MM, NN, [ax, ly]); // grid generation
plot(Sp, ps=“ff_mesh.eps”); // plot and save mesh
//***************************************************************************************************
// ------ Identify Gridpoint Coordinates --------
// Vh xx=x, yy=y;
real[int] xx(Sp.nv);
real[int] yy(Sp.nv);
for(int i = 0; i <Sp.nv; i++){
xx[i] = Sp(i).x;
yy[i] = Sp(i).y;
}
/*
------ Macros --------
/
real sqrt2 = sqrt(2.);
//macro Strain(ur,uz) [dx(ur), ur/x, dy(uz), 0.5(dy(ur)+dx(uz))] //EOM
macro Strain(ur,uz) [dx(ur), ur/x, dy(uz), (dy(ur)+dx(uz))/sqrt2] //EOM
//***************************************************************************************************
// ------ Problem Definition -------
fespace Uh(Sp,[P1,P1]); // piecewise linear FE
fespace Vh(Sp,[P1]);
Uh [ur,uz]; // variational variables
Uh [vr,vz]; // variational variables
Vh xcoord=x;
// Iterate through all the frequencies
for(int jj=0; jj<ff.n; jj++) { // for all frequencies
real f0 = ff[jj], w0 = 2pif0; cout << f0/1e3 << “kHz” << endl;
solve AxiDisc([ur,uz],[vr,vz,]) = // variational equation!
int2d(Sp)(0.5xrhoxcoordw0^2*[vr,vz]'[ur,uz])
-int2d(Sp)(xStrain(vr,vz)'StiffStrain(ur,uz)
)
-int1D(Uh, 2)(
force*uz
)
- on(4, ur=0) // BC: axis of symmetry
- on(4, uz=0); // BC: axis of symmetry
Vh urreal=real(ur);
Vh uzreal=real(uz);
mesh Sp2 = movemesh (Sp,[x+1000.*urreal, y+1000.*uzreal]);
plot(Sp,Sp2,urreal, wait=true, cmm=“f = " + f0/1000 + " kHz”, fill = true, value=true);
// plot(Sp,Sp2,uzreal, wait=true, cmm=“f = " + f0/1000 + " kHz”, fill = true, value=true);
}