Hello, i have written a code for an FSI analysis for a 3d model an i tried to apply also a wave formula to the fluid in my model. When i run the code i see that i don’t get any errors but i also don’t get all the results that i ask such as the matrix of the displacements, the plot of the displacements, the plot of the forces. Can anyone tell me if there is anything wrong with my code
// Load necessary FreeFem++ libraries
load “msh3”
load “medit”
load “gmsh”
// Parameters
int bottombeam = 2; //label of bottombeam
real E = 21.5;
real sigma = 0.29;
real gravity = -0.05;
real coef = 0.2;
// Time parameters
real T = 1.0; // Total time
real dt = 0.01; // Time step
int nsteps = T/dt; // Number of time steps
real freq = 1.0; // Frequency of the wave
real amp = 1.0; // Amplitude of the wave
// Mesh
border a(t=0, 1) {x=t; y=0; z=0; label=1;};
border b(t=0, 2) {x=1; y=t; z=0; label=1;};
border c(t=1, 3) {x=t; y=2; label=1;};
border d(t=0, 2) {x=3; y=2-t; label=1;};
border e(t=0, 1) {x=3+t; y=0; label=1;};
border f(t=0, 1) {x=4; y=3t; label=1;};
border g(t=1, 0) {x=4t; y=3; label=1;};
border h(t=1, 0) {x=0; y=3*t; label=1;};
mesh th1 = buildmesh(a(10) + b(15) + c(20) + d(15) + e(10) + f(20) + g(30) + h(20));
plot(th1, wait=true);
// Construction 3d mesh
func zmin = 0;
func zmax = 2.;
int MaxLayer = 10;
mesh3 th2 = buildlayers(th1, MaxLayer, zbound=[zmin,zmax]);
medit(“th1”, th2);
// Fespace (Solid)
fespace Vh(th2,[P1,P1,P1]);
Vh [uu1,uu2,uu3], [v1,v2,v3];
// Macro
real sqrt2=sqrt(2.);
macro epsilon(u1,u2,u3) [dx(u1),dy(u2),dz(u3),(dz(u2)+dy(u3))/sqrt2,(dz(u1)+dx(u3))/sqrt2,(dy(u1)+dx(u2))/sqrt2] // EOM
macro div(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM
real mu = E/(2*(1+sigma));
real lambda = Esigma/((1+sigma)(1-2sigma));
cout <<"Lambda = "<< gravity << endl;
cout << "Mu = " << mu << endl;
solve Elasticity([uu1,uu2,uu3],[v1,v2,v3])=
int3d(th2)(
lambdadiv(uu1,uu2,uu3)*div(v1,v2,v3)
+2.mu(epsilon(uu1,uu2,uu3)'*epsilon(v1,v2,v3) )
)
- int2d(th2) (gravity*v3)
- on(1,uu1=0,uu2=0,uu3=0)
;
cout << "Max displacement = " << uu1.linfty << endl;
cout << "Max displacement = " << uu2.linfty << endl;
cout << "Max displacement = " << uu3.linfty << endl;
plot([uu1,uu2,uu3], wait=true);
plot(th2, [uu1,uu2,uu3], wait=1, value=true);
// Movemesh
mesh3 th3 = movemesh(th2, [x, y+uu2, z]);
mesh3 th4 = movemesh3(th2, transfo=[x+uu1, y+uu2, z+uu3]);
// Plot
plot(th3, wait=true, value=true);
// Mesh (fluid)
border k(t=0, 4){x=t; y=-1; label= 1;}
border l(t=-1, 1){x=4; y=t ; label= 2;}
border m(t=0, 4){x=4-t; y=1; label= 3;}
border n(t=1, -1){x=0; y=t; label= 4;}
mesh sh1 = buildmesh(k(10) + l(10) + m(10) + n(10));
plot(sh1, wait=true);
// Construction 3d mesh (Fluid)
func zmin2 = 0;
func zmax2 = 2.;
int MaxLayer2 = 10;
mesh3 sh2 = buildlayers(sh1, MaxLayer2, zbound=[zmin2,zmax2]);
medit(“sh1”, sh2);
// Mesh (total)
mesh3 qh = th2 + sh2;
plot(qh, wait=1);
medit(“qh”, qh);
// Fespace (fluid)
fespace Xh(sh2, [P2, P2, P2]);
Xh [uf1, uf2, uf3];
Xh [bc1, bc2, bc3];
Xh [brhs1, brhs2, brhs3];
Xh [bcx, bcy, bcz];
fespace Mh(sh2, [P1, P1, P1]);
Mh [p, pp, ppp];
Mh [q, qq, qqq];
// Macro (Fluid)
macro Div(uf1, uf2, uf3) (dx(uf1) + dy(uf2) + dz(uf3)) // EOM
varf stokesPressure(p, q) = int3d(sh2)(q*Div(uf1, uf2, uf3)); // EOM
// Define the variational forms for the velocity components and pressure
varf bx (u1, q) = int3d(sh2)(-(dx(uf1)*q));
varf by (u1, q) = int3d(sh2)(-(dy(uf1)*q));
varf bz (u1, q) = int3d(sh2)(-(dz(uf1)*q));
// Problem (Fluid)
varf Lap ([uf1, uf2, uf3],[vuf1, vuf2, vuf3])
= int3d(sh2)(
dx(uf1)*dx(vuf1) + dy(uf1)*dy(vuf1) + dz(uf1)*dz(vuf1)
+ dx(uf1)*dx(vuf2) + dy(uf1)*dy(vuf2) + dz(uf1)*dz(vuf2)
+ dx(uf1)*dx(vuf3) + dy(uf1)*dy(vuf3) + dz(uf1)*dz(vuf3)
+ dx(uf2)*dx(vuf1) + dy(uf2)*dy(vuf1) + dz(uf2)*dz(vuf1)
+ dx(uf2)*dx(vuf2) + dy(uf2)*dy(vuf2) + dz(uf2)*dz(vuf2)
+ dx(uf2)*dx(vuf3) + dy(uf2)*dy(vuf3) + dz(uf2)*dz(vuf3)
+ dx(uf3)*dx(vuf1) + dy(uf3)*dy(vuf1) + dz(uf3)*dz(vuf1)
+ dx(uf3)*dx(vuf2) + dy(uf3)*dy(vuf2) + dz(uf3)*dz(vuf2)
+ dx(uf3)*dx(vuf3) + dy(uf3)*dy(vuf3) + dz(uf3)*dz(vuf3)
)
+ on(2, uf1=1, uf2=1, uf3=1)
+ on(1,3, uf1=0, uf2=0, uf3=0)
;
plot(sh2, [uf1, uf2, uf3], wait=1, value=true);
varf Mass(p,q)=int2d(sh1)(p*q);
matrix A = Lap(Xh, Xh, solver=CG);
matrix M = Mass(Mh, Mh,solver=CG);
matrix Bx = bx(Xh, Mh);
matrix By = by(Xh, Mh);
matrix Bz = bz(Xh, Mh);
func real[int] divup (real[int] & pp){
int verb = verbosity;
verbosity = 1;
brhs1 = Bx’pp; brhs1[] += bc1[] .bcx[];
uf1[] = A^-1brhs1[];
brhs2[] = By’pp; brhs2[] += bc2[] .bcy[];
uf2[] = A^-1brhs2[];
brhs3[] = Bz’pp; brhs3[] += bc3[] .bcz[];
uf3[] = A^-1brhs3[];
ppp[] = M^-1pp;
ppp[] = 1.e-6ppp[];
ppp[] = Bxuf1;
ppp += Byuf2[];
ppp[] += Bzuf3;
verbosity=verb;
return ppp ;
};
p=0;q=0;uf1=0;uf2=0;uf3=0;v1=0;v2=0;v3=0;
// Coupling loop
real time = 0;
for(int step = 0; step < nsteps; ++step){
time += dt;
// Apply time-dependent boundary condition
varf LapWave ([uf1, uf2, uf3], [vuf1, vuf2, vuf3])
= int3d(sh2)(
dx(uf1)*dx(vuf1) + dy(uf1)*dy(vuf1) + dz(uf1)*dz(vuf1)
+ dx(uf2)*dx(vuf2) + dy(uf2)*dy(vuf2) + dz(uf2)*dz(vuf2)
+ dx(uf3)*dx(vuf3) + dy(uf3)*dy(vuf3) + dz(uf3)*dz(vuf3)
)
+ on(2, uf1=amp * sin(2*pi*freq*time), uf2=amp * sin(2*pi*freq*time), uf3=amp * sin(2*pi*freq*time))
+ on(1,3, uf1=0, uf2=0, uf3=0);
matrix AWave = LapWave(Xh, Xh, solver=CG);
real[int] rhsWave = LapWave(0, Xh);
uf1[] = AWave^-1 * rhsWave;
plot(sh2, [uf1, uf2, uf3], wait=1, value=true);
plot(sh2, wait=1, value=true);
// Solve (fluid)
LinearCG(divup, p[], eps=1.e-3, nbiter=50);
divup(p[]);
// Forces
Vh [sigma11, sigma22, sigma33];
Vh [sigma12, sigma13, sigma23];
Vh [u1, u2, u3] = [uu1, uu2, uu3];
sigma11 = lambda*(dx(u1) + dy(u2) + dz(u3)) + 2mudx(u1);
sigma22 = lambda*(dx(u1) + dy(u2) + dz(u3)) + 2mudy(u2);
sigma33 = lambda*(dx(u1) + dy(u2) + dz(u3)) + 2mudz(u3);
sigma12 = mu*(dy(u1) + dx(u2));
sigma13 = mu*(dz(u1) + dx(u3));
sigma23 = mu*(dz(u2) + dy(u3));
// Solve (solid)
solve Elasticity2 ([uu1, uu2, uu3], [v1, v2, v3], init=step)
= int3d(th2)(
lambda*div(v1,v2,v3)*div(uu1,uu2,uu3)
+ 2.mu(epsilon(v1,v2,v3)'*epsilon(uu1,uu2,uu3))
)
- int3d(th2)(
- gravity*v3
)
- gravity*v3
- int3d(th2)(
- coef*(sigma11N.xv1 + sigma22N.yv2 + sigma33N.zv3 + sigma12*(N.yv1+N.xv2) + sigma13*(N.zv1+N.xv3) + sigma23*(N.zv2+N.yv3))
)
- coef*(sigma11N.xv1 + sigma22N.yv2 + sigma33N.zv3 + sigma12*(N.yv1+N.xv2) + sigma13*(N.zv1+N.xv3) + sigma23*(N.zv2+N.yv3))
- on(1,uu1=0,uu2=0,uu3=0)
;
// Plot
plot([uu1, uu2, uu3], value=true, cmm=“u”, wait=0);
plot([uf1, uf2, uf3], value=true, cmm=“uf”, wait=0);
// Plot stress components
plot(sigma11, value=true, cmm=“sigma11”, wait=0);
plot(sigma22, value=true, cmm=“sigma22”, wait=0);
plot(sigma33, value=true, cmm=“sigma33”, wait=0);
plot(sigma12, value=true, cmm=“sigma12”, wait=0);
plot(sigma13, value=true, cmm=“sigma13”, wait=0);
plot(sigma23, value=true, cmm=“sigma23”, wait=0);
// Error
real err = sqrt(int3d(th2)((uu1-uf1)^2 + (uu2-uf2)^2 + (uu3-uf3)^2));
cout << "Erreur L2 = " << err << endl;
// Movemesh
th2 = movemesh(th2, [x+0.2uu1, y+0.2uu2, z+0.2*uu3]);
plot(th2, wait=true);
};