0
$\begingroup$

I'm testing a Finite Elements code. To do this, i have created a simple one element quadrilateral mesh to verify the solution:

enter image description here

This is a plane stress elasticity problem, and the local stiffness is computed through the following weak formulation:

enter image description here

The problem data is the following:

enter image description here

the resultant stiffness matrix and load vector:

enter image description here

enter image description here

After applying a penalty boundary conditions (multiply key boundary points in the matrix by a very large number) the final stiffness matrix is:

enter image description here

and after solving the linear system above, the solution is:

enter image description here

Now, i'm trying to solve the same problem in a rotated mesh (45 degrees counter clockwise), and i'm facing some problems to impose the dirichlet boundary conditions, i.e. to impose the null displacement in the specified nodes.

Here are the new nodes coordnates:

enter image description here

the new mesh is:

enter image description here

the rotated stiffness matrix is the following:

enter image description here

and the rotated load vector is:

enter image description here

The solution must be the same in both meshes, but i'm facing a problem to impose the new boundary conditions. How can i impose the natural boundary conditions in the rotated mesh?

Any help will be appreciated. Thank you.

$\endgroup$
  • $\begingroup$ You need to specify the equations you want to solve, along which kind of method you use to discretize the equations onto an element. Otherwise, we do not know what exactly you want to do. $\endgroup$ – Wolfgang Bangerth Jan 17 '17 at 1:00
  • $\begingroup$ I have edited the problem including the formulation to compute the local stiffness matrix K. The method used to reach the weak formulation is the Galerkin method. $\endgroup$ – Diogo Jan 17 '17 at 1:24
  • $\begingroup$ the boundary conditions have changed when the mesh is rotated. $\endgroup$ – mronlinetutor Jan 22 '17 at 19:32
1
$\begingroup$

The large-number approach you are using to impose constraints is generally referred to as a penalty method. One upside of penalty methods is that they are generally easy to apply. A downside is that adding large numbers to the equations makes them ill-conditioned and can cause errors in the solution.

One way to look at penalty methods is to write the total potential energy of your system

$$ E_{TP} = \frac{1}{2}U^TKU - U^TF $$

where $U$ is the vector of nodal displacments, $K$ is the global stiffness matrix, and $F$ is the load vector. Constraints on the values of $U$ can be prescribed by including an additional term in the potential energy

$$ E_{TP} = \frac{1}{2}U^TKU - U^TF + \frac{r}{2}\sum_{i=1}^{n_c} g_i(U)^2 $$

where the $g_i(U)$ are constaint equations that are functions of $U$, $r$ is a large number, and $n_c$ is the number of constraints. The equations, $g_i$, equal zero when the constraint is satisfied. In your un-rotated case, the $g_i$ are simply the components of $U_i$ you want to constrain.

In the rotated case, you want to impose constraints on the displacements in the rotated coordinate system. At a particular node $j$ these can be written as

$$ \left\{\begin{array}{c}u^\prime_j \\ v^\prime_j\end{array}\right\} = \left[\begin{array}{cc}cos\theta & sin\theta \\ -sin\theta & cos\theta\end{array}\right] \left\{\begin{array}{c}u_j \\ v_j\end{array}\right\} $$

So a typical $g_i$ might be $cos\theta u_j + sin\theta v_j$ if you wanted to constrain $u^\prime_j$ to zero. By minimizing the potential energy with respect to the $U$ components, just like you did to find the stiffness matrix and load vector, there will also be a set of $2 \times 2$ constraint matrices, one for each displacement constraint you want to apply. These are added to the global stiffness matrix just like the large diagonal term you added in the un-rotated problem.

You can find more information on applying constraints in FEM in these sets of notes from Carlos Felippa

MultiFreedom Constraints I

MultiFreedom Constraints II

He refers to these types of constraints as "multifreedom" because they are a linear combination of more than a single degree-of-freedom in the FE model. He also discusses the alternate, Lagrange multiplier approach which is a better approach, numerically, to apply these types of constraints.

$\endgroup$
  • $\begingroup$ Sorry, i'm not an specialist in finite elemts, but anyway thank you for your answer, it is of great help. Even whith your help i still cant solve the problem. If it is not a problem to you, could you show the final numerical rotated stiffness matrix whith its boundary conditions? I belive it would be easier to me to understand how to solve. $\endgroup$ – Diogo Jan 18 '17 at 11:20
  • $\begingroup$ I don't see how having numerical values for this specific problem would be useful. However, I did add links to sets of notes that explain FE constraints in more detail. $\endgroup$ – Bill Greene Jan 18 '17 at 17:20
0
$\begingroup$

if u impose a boundary condition u1=v1=0, v2=0, u4=0, the single rectangular mesh shouldn't rotate at all because node one is locked into place, node two can stretch horizontally while node 4 can stretch vertically leaving node three free to mode. applying a stretching force of 6000lbs in the horizontal direction at node two and at node three will deform the elements from a rectangular shape to some unknown shape but not rotate the rectangular mesh.

when the rectangular mesh is rotated, the u and v boundary conditions have changed even though you have rotated the input forces. the stiffness matrix should remain the same only the boundary conditions have changed.

[K]=[J]*integral [B]^t[D][B]dV where [F]=[K][x]. rotating [x] whould change the integral by a factor equals to the jacobian matrix but not the stiffness matrix. [F]=[K][x], [e]=[Q][u], [u]=[N][d], [e]=[B][d],where Q is the differential operator.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.