我遵循本教程,以尝试和做一个model.msh
的有限元分析,我必须看看它会如何变形给定不同的外力在不同的地方。
他们把弱形式定义为
a(u,v) = ∫( ε(v) ⊙ (σ∘ε(u)) )*dΩ
l(v) = 0
,它们表示"线性形式简单为l(v) = 0,因为在这个例子中没有外力。">
如前所述,我想分析不同外力会对我的模型造成的不同变形,但我似乎找不到任何这样的例子。谁能帮我定义一下大于0的外力的线性形式?
谢谢。
也许这对你有帮助。这是仓促写的,所以如果你遇到拼写错误或其他美观问题,请不要介意:)
# define where the output shall go
output_Path ="Output/3_Gridap/1_Lin_FEA/FE_8"
mkpath(output_Path)
output_Name ="pde_6"
using Gridap
# please load the model that is shown in: https://gridap.github.io/Tutorials/dev/pages/t001_poisson/
model = DiscreteModelFromFile("Code/Meshes/Data/possion.json")
# just in case you want to see the model using paraview
writevtk(model,"$(output_Path)/md")
order = 1
reffe = ReferenceFE(lagrangian,VectorValue{3,Float64},order)
V0 = TestFESpace(model,reffe;
conformity=:H1,
# to see which elements belongs to "bottom" open the model which is saved through "writevtk(model,"$(output_Path)/md")"
dirichlet_tags=["bottom"],
# activate/deactivate the boundary conditions
dirichlet_masks=[
(true, true, true), # clamp the bottom
])
# define displacement
clamping(x) = VectorValue(0.0,0.0,0.0)
U = TrialFESpace(V0,[clamping])
const E = 7e+7
const ν = 0.33
const λ = (E*ν)/((1+ν)*(1-2*ν))
const μ = E/(2*(1+ν))
σ(ε) = λ*tr(ε)*one(ε) + 2*μ*ε
degree = 2*order
Ω = Triangulation(model)
dΩ = Measure(Ω,degree)
# Neumann boundary conditions
# here we define the surface on which we want an external force
Γ_Tr = BoundaryTriangulation(model,tags=["triangle"])
dΓ_Tr = Measure(Γ_Tr,degree)
# a force shall be applied on the y-direction
f_Tr(x) = VectorValue(0.0, 1e+6, 0.0)
# mass forces due to gravity, the value is set quite high, such that an impact can be seen
mass_Forces(x) = VectorValue(0.0, -1e+7, 0.0)
# Weak form
a(u,v) = ∫( ε(v) ⊙ (σ∘ε(u)) )*dΩ
l(v) = ∫( v ⋅ mass_Forces )* dΩ + ∫( v ⋅ f_Tr )* dΓ_Tr
op = AffineFEOperator(a,l,U,V0)
uh = solve(op)
writevtk(Ω,"$(output_Path)/$(output_Name)",
cellfields=[
"uh" => uh,
"epsi" => ε(uh),
"sigma" => σ∘ε(uh)])