# Fluid-Structure Analysis 3d model

**URL:** <https://community.freefem.org/t/fluid-structure-analysis-3d-model/3265>\
**Category:** General Discussion\
**Created:** [May 23, 2024, 10:20am UTC](https://community.freefem.org/t/fluid-structure-analysis-3d-model/3265 "2024-05-23T10:20:28Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![zoumpou](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/zoumpou/32/2214_2.png) [@zoumpou](https://community.freefem.org/u/zoumpou)\
**Post date:** [May 23, 2024, 10:20am UTC](https://community.freefem.org/t/fluid-structure-analysis-3d-model/3265/1 "2024-05-23T10:20:28Z")

</div>

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=3_t; label=1;};  
border g(t=1, 0) {x=4_t; 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 = E_sigma/((1+sigma)_(1-2_sigma));  
cout \<\<"Lambda = "\<\< gravity \<\< endl;  
cout \<\< "Mu = " \<\< mu \<\< endl;  
solve Elasticity([uu1,uu2,uu3],[v1,v2,v3])=  
int3d(th2)(  
lambda_div(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[] = Bx_uf1;  
ppp += By_uf2[];  
ppp[] += Bz_uf3;  
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)) + 2_mu_dx(u1);  
sigma22 = lambda\*(dx(u1) + dy(u2) + dz(u3)) + 2_mu_dy(u2);  
sigma33 = lambda\*(dx(u1) + dy(u2) + dz(u3)) + 2_mu_dz(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  
)

- int3d(th2)(
  - coef\*(sigma11_N.x_v1 + sigma22_N.y_v2 + sigma33_N.z_v3 + sigma12\*(N.y_v1+N.x_v2) + sigma13\*(N.z_v1+N.x_v3) + sigma23\*(N.z_v2+N.y_v3))  
)

- 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.2_uu1, y+0.2_uu2, z+0.2\*uu3]);  
plot(th2, wait=true);  
};