Problem
So I have a simple rectangular domain with x = 1 m and y = 0.5 m. What I want to compute is a steady state if pressure P is applied to the left and right vertical side.
Attempt at solution
load "medit"
// Generate mesh.
real x0 = 0;
real x1 = 1;
real y0 = 0;
real y1 = 0.5;
int n = 20;
real m = 10;
mesh Th = square(n, m, [x0+(x1-x0)*x, y0+(y1-y0)*y]);
plot(Th, wait=1, cmm="Mesh");
// **********
// Variables.
// **********
real sqrt2 = sqrt(2.);
real totalTime = 150.; // Total simulation time.
real dt = 5.0; // Time step.
real p = 1e8;
// Material properties.
real youngModul = 72.1e9; // Elasticity (Young) modulus.
real poisson = 0.33; // Poisson ratio.
// Lame parameters
real mu = youngModul / 2. / (1. + poisson);
real lam = youngModul * poisson / (1. - 2. * poisson) / (1. + poisson);
real G = youngModul / 2. / (1. + poisson);
// *********************
// Finite element space.
// *********************
// The finite element space defined over mesh
fespace Vh(Th, P1);
fespace Nh(Th, [P1, P1]);
fespace Zh(Th, [P1, P1, P1]);
// FEM variables.
Vh v;
Nh [u1,u2], [v1,v2];
Zh [Sxx, Syy, Sxy];
// *******
// Macros.
// *******
macro stress [Sxx, Syy, Sxy] // stress tensor
macro grad(u) [dx(u), dy(u)] //EOM
macro div(u1, u2) (dx(u1) + dy(u2)) // EOM
macro epsilon(u1, u2) [dx(u1), dy(u2), (dy(u1) + dx(u2)) / sqrt2] // EOM
macro e(u1, u2) [dx(u1), dy(u2), dz(u3), dy(u1) + dx(u2)] // EOM
macro matE [[(lam + 2. * mu), lam, 0.],
[lam, (lam + 2. * mu), 0.],
[0., 0., mu]] // EOM plane strain stress-strain matrix
So far, things should be ok. Now I am having troubles defining the problem
// ***********************
// Variational definition.
// ***********************
varf elas([u1,u2],[v1,v2]) =
int2d(Th)(
lambda*div(u1,u2)*div(v1,v2)
+2.*mu*( epsilon(u1,u2)'*epsilon(v1,v2) )
)
;
varf l([u1,u2],[v1,v2]) = ;// on(2,u1=0,u2=0) what?
matrix A = elas(Vh, Vh); //stiffness matrix
because I can’t figure out how to set the boundary condition. If the boundary condition was deflection then it is obvious, but now I do not know hot to change the problem definition so that the label 2 would have the 1e8 pressure in the x direction and label 4 the same amount in the opposite direction.