Basically as the title says Im learning FreeFem for a university project and I want to better understand what the program is actually doing. I have the following program modified from https://modules.freefem.org/modules/elasticity/:
load "msh3" load "gmsh" //* mesh int Clamped = 27; //Beam fixed mesh3 Th = readmesh3("rails"); plot(Th, wait=true); //* Some parameters defined here... real E = 210e9; //* Young's Modulus real Nu = 0.27; //* poisson's ratio real Rho = 7800; //* density real Gravity = -9.81; //* gravitational acceleration //* fespace func Pk = P1; fespace Uh(Th, [Pk, Pk, Pk]); Uh [ux, uy, uz]; Uh [vx, vy, vz]; Uh [uxp, uyp, uzp]; Uh [uxpp, uypp, uzpp]; //* Macro real sqrt2 = sqrt(2.); macro Epsilon(ux, uy, uz) [dx(ux), dy(uy), dz(uz), (dz(uy)+dy(uz))/sqrt2, (dz(ux)+dx(uz))/sqrt2, (dy(ux)+dx(uy))/sqrt2] // macro Divergence(ux, uy, uz) (dx(ux) + dy(uy) + dz(uz)) // //* Problem real Mu = E/(2.*(1.+Nu)); real Lambda = E*Nu/((1.+Nu)*(1.-2.*Nu)); varf vElasticity ([ux, uy, uz], [vx, vy, vz]) = int3d(Th)( Lambda * Divergence(vx, vy, vz) * Divergence(ux, uy, uz) + 2. * Mu * ( Epsilon(vx, vy, vz)' * Epsilon(ux, uy, uz) ) ) + int3d(Th)( Rho * Gravity * vz ) - int2d(Th, 30)( 30000*vy ) + on(Clamped, ux=0, uy=0, uz=0) + on(26, ux=0, uy=0, uz=0) ; matrix Elasticity = vElasticity(Uh, Uh, solver=sparsesolver); real[int] ElasticityBoundary = vElasticity(0, Uh); ux = Elasticity^-1 * ElasticityBoundary; //*Movemesh real dmax = ux.max; cout << "max displacement = " << dmax << endl; real coef = 50; //Th = movemesh(Th, [x+ux, y+uy, z+uz]); //[ux, uy, uz] = [ux, uy, uz]; mesh3 Thm = movemesh3(Th, transfo=[x+ux*coef, y+uy*coef, z+uz*coef]); Thm = change(Thm); //*Plot //plot([ux, uy, uz], value=true); plot(Th, Thm, value=true, wait=true, cmm="coef amplification = "+coef);
A few of the points Im confused about are:
under fespace, what does
Uh [uxp, uyp, uzp]and
Uh [uxpp, uypp, uzpp]do/represent? These aren’t used elsewhere so does not sure what it does.
when solving the matrix vector equation Ax=b, why do the results get written to
ux? I have applied the forces in x,y and z direction but it doesn’t seem to matter? so why the x in
When applying external forces I am using
int2d(Th, No.). What does
int2drepresent? I also found
int0dfor example be a point force? And the extra number I have to pass, does that correspond to a line, point or surface on my mesh? (Im preparing in gmsh fyi).
The program I have can display the displacements and not the stresses. I am actually mostly interested in these stresses. What modifications would I need to calculate the stresses in the structure?
Sorry if my code is a bit bad, Im very new to this software and am trying to understand it better. Any help is greatly appreciated!
p.s. This is my result so far: