In [1]:
import getfem as gf
import numpy as np
In [2]:
#
# Physical parameters
#
epsilon = 1. # Thickness of the plate (cm)
E = 21E6 # Young Modulus (N/cm^2)
nu = 0.3 # Poisson ratio
clambda = E*nu/((1+nu)*(1-2*nu)) # First Lame coefficient (N/cm^2)
cmu = E/(2*(1+nu)) # Second Lame coefficient (N/cm^2)
clambdastar = 2*clambda*cmu/(clambda+2*cmu) # Lame coefficient for Plane stress (N/cm^2)
F = 100E2 # Force density at the right boundary (N/cm^2)
kappa = 4. # Thermal conductivity (W/(cm K))
D = 10. # Heat transfert coefficient (W/(K cm^2))
air_temp = 20. # Temperature of the air in oC.
alpha_th = 16.6E-6 # Thermal expansion coefficient (/K).
T0 = 20. # Reference temperature in oC.
rho_0 = 1.754E-8 # Resistance temperature coefficient at T0 = 20oC
alpha = 0.0039 # Second resistance temperature coefficient.
#
# Numerical parameters
#
h = 2. # Approximate mesh size
elements_degree = 2 # Degree of the finite element methods
export_mesh = True # Draw the mesh after mesh generation or not
solve_in_two_steps = True # Solve the elasticity problem separately or not
In [3]:
mo1 = gf.MesherObject('rectangle', [0., 0.], [100., 25.])
mo2 = gf.MesherObject('ball', [25., 12.5], 8.)
mo3 = gf.MesherObject('ball', [50., 12.5], 8.)
mo4 = gf.MesherObject('ball', [75., 12.5], 8.)
mo5 = gf.MesherObject('union', mo2, mo3, mo4)
mo = gf.MesherObject('set minus', mo1, mo5)
In [4]:
mesh = gf.Mesh('generate', mo, h, 2)
In [5]:
fb1 = mesh.outer_faces_in_box([1., 1.], [99., 24.])
fb2 = mesh.outer_faces_with_direction([ 1., 0.], 0.01)
fb3 = mesh.outer_faces_with_direction([-1., 0.], 0.01)
fb4 = mesh.outer_faces_with_direction([0., 1.], 0.01)
fb5 = mesh.outer_faces_with_direction([0., -1.], 0.01)
RIGHT_BOUND=1; LEFT_BOUND=2; TOP_BOUND=3; BOTTOM_BOUND=4; HOLE_BOUND=5;
mesh.set_region( RIGHT_BOUND, fb2)
mesh.set_region( LEFT_BOUND, fb3)
mesh.set_region( TOP_BOUND, fb4)
mesh.set_region(BOTTOM_BOUND, fb5)
mesh.set_region( HOLE_BOUND, fb1)
mesh.region_subtract( RIGHT_BOUND, HOLE_BOUND)
mesh.region_subtract( LEFT_BOUND, HOLE_BOUND)
mesh.region_subtract( TOP_BOUND, HOLE_BOUND)
mesh.region_subtract(BOTTOM_BOUND, HOLE_BOUND)
In [6]:
mesh.export_to_vtk('mesh.vtk');
In [7]:
%%writefile plot.py
from paraview.simple import *
paraview.simple._DisableFirstRenderCameraReset()
mesh_vtk = LegacyVTKReader( FileNames=['mesh.vtk'] )
RenderView1 = GetRenderView()
RenderView1.CameraFocalPoint = [50.0, 12.5, 0.0]
RenderView1.CameraPosition = [50.0, 12.5, 10000.0]
RenderView1.InteractionMode = '2D'
RenderView1.CenterOfRotation = [50.0, 12.5, 0.0]
DataRepresentation1 = Show()
DataRepresentation1.ScaleFactor = 10.0
DataRepresentation1.ScalarOpacityUnitDistance = 10.363332599018701
DataRepresentation1.EdgeColor = [0.0, 0.0, 0.5000076295109483]
RenderView1.CameraPosition = [50.0, 12.5, 199.13071041509224]
RenderView1.CenterAxesVisibility = 0
RenderView1.CameraClippingRange = [197.1394033109413, 202.11767107131863]
RenderView1.CameraParallelScale = 51.53882032022076
DataRepresentation1.Representation = 'Wireframe'
WriteImage('mesh.png')
Render()
In [8]:
!python plot.py
In [9]:
from IPython.core.display import Image
Image('mesh.png')
Out[9]:
In [10]:
mfu = gf.MeshFem(mesh, 2) # Finite element for the elastic displacement
mfu.set_classical_fem(elements_degree)
mft = gf.MeshFem(mesh, 1) # Finite element for temperature and electrical field
mft.set_classical_fem(elements_degree)
mfvm = gf.MeshFem(mesh, 1) # Finite element for Von Mises stress interpolation
mfvm.set_classical_discontinuous_fem(elements_degree)
mim = gf.MeshIm(mesh, pow(elements_degree,2)) # Integration method
In [11]:
md=gf.Model('real');
md.add_fem_variable('u', mfu) # Displacement of the structure
md.add_fem_variable('theta', mft) # Temperature
md.add_fem_variable('V', mft) # Electric potential
In [12]:
md.add_initialized_data('cmu', [cmu])
md.add_initialized_data('clambdastar', [clambdastar])
md.add_isotropic_linearized_elasticity_brick(mim, 'u', 'clambdastar', 'cmu')
Out[12]:
In [13]:
md.add_Dirichlet_condition_with_multipliers(mim, 'u', elements_degree-1, LEFT_BOUND)
md.add_initialized_data('Fdata', [F*epsilon, 0])
md.add_source_term_brick(mim, 'u', 'Fdata', RIGHT_BOUND)
Out[13]:
In [14]:
sigmaeps = '(eps/(rho_0*(1+alpha*(theta-T0))))'
md.add_initialized_data('eps', [epsilon])
md.add_initialized_data('rho_0', [rho_0])
md.add_initialized_data('alpha', [alpha])
md.add_initialized_data('T0', [T0])
md.add_nonlinear_generic_assembly_brick(mim, sigmaeps+'*(Grad_V.Grad_Test_V)')
md.add_Dirichlet_condition_with_multipliers(mim, 'V', elements_degree-1, RIGHT_BOUND)
md.add_initialized_data('DdataV', [0.1])
md.add_Dirichlet_condition_with_multipliers(mim, 'V', elements_degree-1, LEFT_BOUND, 'DdataV')
Out[14]:
In [15]:
md.add_initialized_data('kaeps', [kappa*epsilon])
md.add_generic_elliptic_brick(mim, 'theta', 'kaeps')
md.add_initialized_data('D2', [D*2])
md.add_initialized_data('D2airt', [air_temp*D*2])
md.add_mass_brick(mim, 'theta', 'D2')
md.add_source_term_brick(mim, 'theta', 'D2airt')
md.add_initialized_data('Deps', [D/epsilon])
md.add_initialized_data('Depsairt', [air_temp*D/epsilon])
md.add_Fourier_Robin_brick(mim, 'theta', 'Deps', TOP_BOUND)
md.add_source_term_brick(mim, 'theta', 'Depsairt', TOP_BOUND)
md.add_Fourier_Robin_brick(mim, 'theta', 'Deps', BOTTOM_BOUND)
md.add_source_term_brick(mim, 'theta', 'Depsairt', BOTTOM_BOUND)
Out[15]:
In [16]:
md.add_nonlinear_generic_assembly_brick(mim, '-'+sigmaeps+'*Norm_sqr(Grad_V)*Test_theta')
Out[16]:
In [17]:
md.add_initialized_data('beta', [alpha_th*E/(1-2*nu)])
md.add_linear_generic_assembly_brick(mim, 'beta*(T0-theta)*Trace(Grad_Test_u)')
Out[17]:
In [18]:
if (solve_in_two_steps):
md.disable_variable('u')
print 'First problem with', md.nbdof(), ' dofs'
md.solve('max_res', 1E-9, 'max_iter', 100, 'noisy')
md.enable_variable('u')
md.disable_variable('theta')
md.disable_variable('V')
print 'Second problem with ', md.nbdof(), ' dofs'
md.solve('max_res', 1E-9, 'max_iter', 100, 'noisy')
else:
print 'Global problem with ', md.nbdof(), ' dofs'
md.solve('max_res', 1E-9, 'max_iter', 100, 'noisy')
In [19]:
U = md.variable('u')
V = md.variable('V')
THETA = md.variable('theta')
VM = md.compute_isotropic_linearized_Von_Mises_or_Tresca('u', 'clambdastar', 'cmu', mfvm)
CO = np.reshape(md.interpolation('-'+sigmaeps+'*Grad_V', mfvm), (2, mfvm.nbdof()), 'F')
mfvm.export_to_vtk('displacement_with_von_mises.vtk', mfvm, VM, 'Von Mises Stresses', mfu, U, 'Displacements')
mft.export_to_vtk('temperature.vtk', mft, THETA, 'Temperature')
mft.export_to_vtk('electric_potential.vtk', mft, V, 'Electric potential')
In [20]:
%%writefile plot.py
from paraview.simple import *
paraview.simple._DisableFirstRenderCameraReset()
displacement_with_von_mises_vtk = LegacyVTKReader( FileNames=['displacement_with_von_mises.vtk'] )
RenderView1 = GetRenderView()
RenderView1.CameraFocalPoint = [50.0, 12.5, 0.0]
RenderView1.CameraPosition = [50.0, 12.5, 10000.0]
RenderView1.InteractionMode = '2D'
RenderView1.CenterOfRotation = [50.0, 12.5, 0.0]
DataRepresentation1 = Show()
DataRepresentation1.EdgeColor = [0.0, 0.0, 0.5000076295109483]
DataRepresentation1.SelectionPointFieldDataArrayName = 'Von_Mises_Stresses'
DataRepresentation1.ColorArrayName = ('POINT_DATA', 'Von_Mises_Stresses')
DataRepresentation1.ScalarOpacityUnitDistance = 10.363332599018701
DataRepresentation1.ScaleFactor = 10.0
a1_Von_Mises_Stresses_PVLookupTable = GetLookupTableForArray( "Von_Mises_Stresses", 1, RGBPoints=[16.11453628540039, 0.23, 0.299, 0.754, 24018.4283618927, 0.865, 0.865, 0.865, 48020.7421875, 0.706, 0.016, 0.15], VectorMode='Magnitude', NanColor=[0.25, 0.0, 0.0], ColorSpace='Diverging', ScalarRangeInitialized=1.0 )
a1_Von_Mises_Stresses_PiecewiseFunction = CreatePiecewiseFunction( Points=[16.11453628540039, 0.0, 0.5, 0.0, 48020.7421875, 1.0, 0.5, 0.0] )
DataRepresentation1.ScalarOpacityFunction = a1_Von_Mises_Stresses_PiecewiseFunction
DataRepresentation1.LookupTable = a1_Von_Mises_Stresses_PVLookupTable
a1_Von_Mises_Stresses_PVLookupTable.ScalarOpacityFunction = a1_Von_Mises_Stresses_PiecewiseFunction
RenderView1.CameraPosition = [50.0, 12.5, 199.13071041509224]
RenderView1.CenterAxesVisibility = 0
RenderView1.CameraClippingRange = [197.1394033109413, 202.11767107131863]
RenderView1.CameraParallelScale = 51.53882032022076
WriteImage('von_mises.png')
Render()
In [21]:
!python plot.py
In [22]:
from IPython.core.display import Image
Image('von_mises.png')
Out[22]:
In [ ]: