Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

Boundary conditions for "supported" edge

More
4 years 8 months ago #1834 by jfk68
We try to calculate stress and deformation of rectangular plates under pressure load.

For fixed edges (no displacement, no rotation) we use dirichlet boundary conditions. No problem.

Assuming that we have a plate with some "supported" edges, i.e. elements are restricted in displacement only in a certain direction, e.g. z-axis while rotation and displacement in x/y is allowed.

How could that be specified in .py code. My assumption is that something of the dof of the relevant elements has to be modified, the question is what and how?

The "support" constraints should be applied to the "supported" bc in the model below
Code:
import ngsolve import netgen from netgen.csg import * from netgen.meshing import * from ngsolve import * import ngsolve.internal a = 300.0 b = 100.0 t = 8.0 pressure = 0.01 cube = OrthoBrick( Pnt(-a/2 - 5, -b/2 - 5, -t), Pnt(a/2 + 5, b/2 + 5, 5) ) cube *= Plane(Pnt(0,0,0),Vec(0,0,1)).bc("force") cube *= Plane(Pnt(-a/2,0,0),Vec(-1,0,0)).bc("supported") cube *= Plane(Pnt(a/2,0,0),Vec(1,0,0)).bc("supported") cube *= Plane(Pnt(0,-b/2,0),Vec(0,-1,0)).bc("wall") cube *= Plane(Pnt(0,b/2,0),Vec(0,1,0)).bc("wall") geo = CSGeometry() geo.Add (cube) mesh = Mesh(geo.GenerateMesh(maxh=3)) # What should be done to supported bc ? fes = H1(mesh, order=2, dirichlet="wall", dim=mesh.dim) u,v = fes.TrialFunction(), fes.TestFunction() # linear strain tensor def epsilon(u): return 0.5 * (u.Deriv().trans + u.Deriv())# material parameters E = 1200 nu = 0.35 # lame parameters mu = 1/2*(E / (1+nu)) lam = E * nu / ((1+nu)*(1-2*nu)) a = BilinearForm(fes) a += SymbolicBFI( 2*mu*InnerProduct(epsilon(u),epsilon(v)) + lam*Trace(u.Deriv())*Trace(v.Deriv())) a.Assemble() force = CoefficientFunction( (0,0,-pressure) ) f = LinearForm(fes) f += SymbolicLFI( force*v, definedon=mesh.Boundaries("force")) f.Assemble() u = GridFunction(fes,name="displacement") u.vec.data = a.mat.Inverse(fes.FreeDofs())*f.vec Draw(u) # linear stress tensor def stress(u): return 2*mu*epsilon(u) + lam*Trace(epsilon(u))*Id(mesh.dim) Draw(stress(u),mesh,"stress")
More
4 years 7 months ago #1838 by mneunteufel
Hi jfk68,

for supported boundaries you can "split" the H1 space into its (Cartesian) coordinates and prescribe different Dirichlet boundary conditions
Code:
fesx = H1(mesh, order=2, dirichlet="wall") fesy = H1(mesh, order=2, dirichlet="wall") fesz = H1(mesh, order=2, dirichlet="wall|supported") fes = FESpace( [fesx,fesy,fesz] ) (ux,uy,uz),(vx,vy,vz) = fes.TnT() u = CoefficientFunction( (ux,uy,uz) ) v = CoefficientFunction( (vx,vy,vz) ) gradu = CoefficientFunction( (Grad(ux), Grad(uy), Grad(uz)), dims=(mesh.dim,mesh.dim) ) gradv = CoefficientFunction( (Grad(vx), Grad(vy), Grad(vz)), dims=(mesh.dim,mesh.dim) )

Best regards,
Michael
The following user(s) said Thank You: jfk68
More
4 years 7 months ago #1839 by jfk68
@Michael: big thank you, that helps a lot.

Is there some place in the documentation where I could have found this hidden treasure features on my own?


CU

Juergen
More
4 years 7 months ago #1849 by mneunteufel
I did not find an example in the documentation involving different boundary conditions for different components of a VectorH1.

Best regards,
Michael
Time to create page: 0.125 seconds