StartMula
What to doApa perlu dibuat Follow the plan: before class, in class, after class. Each item links to where it is done and ticks itself when it is. Continue takes you to the first thing not yet done.Ikut pelan: sebelum kelas, dalam kelas, selepas kelas. Setiap item dipautkan ke tempat ia dibuat dan ditanda sendiri apabila selesai. Sambung membawa anda ke perkara pertama yang belum selesai.
Counts toward progress:Dikira dalam kemajuan: the bar averages six parts equally: Notes, Slides, Examples, Exercises, Explorer and the can-do checklist. Each tab’s badge shows its own count.bar ini purata enam bahagian secara sama rata: Nota, Slides, Contoh, Latihan, Explorer dan senarai semak boleh-buat. Lencana pada setiap tab menunjukkan kiraannya sendiri.
Your planPelan anda
- A1 [B] at a point · A2 s₁, s₂, s₃ of the cantilever elementA1 [B] pada satu titik · A2 s₁, s₂, s₃ bagi element cantilever
- A3 the global load vector · A4 assembly entries · A5 reading the solutionA3 global load vector · A4 entri assembly · A5 membaca penyelesaian
- A6 the coarse-mesh error and the mesh for 1%A6 ralat mesh kasar dan mesh untuk 1%
- QR Explorer: Cantilever Explorer, task 1Explorer QR: Cantilever Explorer, tugasan 1
○ to do · ✓ done · – done in class, not tracked. Click an item to go to it.○ belum · ✓ selesai · – dibuat dalam kelas, tidak dijejak. Klik item untuk pergi kepadanya.
What this part coversApa yang bahagian ini liputi
Plane Stress Part 1 gave the PDE in u and v, the boundary conditions and the 2 × 8 matrix [N]; this part turns them into numbers, exactly as Heat Transfer Part 2 did for temperature. Galerkin and integration by parts give [k]{û}ᵀ = {q} + {b}; the strain–displacement matrix [B] = [∂]ᵀ[N] puts [k] in the form t∫[B]ᵀ[E][B] dA, whose 64 entries reduce to twelve closed-form values for a rectangle. Edge tractions and self-weight become nodal loads, the elements are assembled through the dof map, the fixed edge is struck, and the text’s cantilever is solved and judged against beam theory and a refined mesh.Plane Stress Part 1 memberi PDE dalam u dan v, boundary condition dan matrix 2 × 8 [N]; bahagian ini menukarkannya kepada nombor, sama seperti Heat Transfer Part 2 lakukan untuk temperature. Galerkin dan integration by parts memberi [k]{û}ᵀ = {q} + {b}; strain–displacement matrix [B] = [∂]ᵀ[N] meletakkan [k] dalam bentuk t∫[B]ᵀ[E][B] dA, dan 64 entrinya menjadi dua belas nilai closed-form bagi segi empat tepat. Traction tepi dan self-weight menjadi nodal load, element di-assemble melalui dof map, tepi yang fixed dibuang, dan cantilever dalam teks diselesaikan serta dinilai berbanding beam theory dan mesh yang di-refine.
After this part you should be able toSelepas bahagian ini, anda sepatutnya boleh
- Apply Galerkin and integration by parts to the plane-stress PDE, identify [k], {q} and {b}, and write [B] = [∂]ᵀ[N] and [k] = t∫[B]ᵀ[E][B] dA.Menggunakan Galerkin dan integration by parts pada PDE plane stress, mengenal pasti [k], {q} dan {b}, dan menulis [B] = [∂]ᵀ[N] serta [k] = t∫[B]ᵀ[E][B] dA.
- Evaluate the 8 × 8 stiffness matrix of a four-node rectangular element from the closed-form entries s1…s12, and check it.Menilai stiffness matrix 8 × 8 bagi four-node rectangular element daripada entri closed-form s1…s12, dan menyemaknya.
- Convert edge tractions and body forces into equivalent nodal loads, assemble the global system with two dofs per node and impose displacement conditions.Menukar traction tepi dan body force kepada nodal load yang setara, meng-assemble sistem global dengan dua dof setiap node dan mengenakan displacement condition.
- Solve the 2 × 2 cantilever, check its reactions, compare it with beam theory and a refined mesh, and explain why the coarse mesh is too stiff.Menyelesaikan cantilever 2 × 2, menyemak reaction-nya, membandingkannya dengan beam theory dan mesh yang di-refine, dan menerangkan mengapa mesh kasar terlalu stiff.
Text: Mohd Yassin et al., Two-Dimensional Finite Element Formulations for Heat, Solids and Fluids (with Matlab), MSNM 2020, §3.2.2–§3.4.Teks: Mohd Yassin et al., Two-Dimensional Finite Element Formulations for Heat, Solids and Fluids (with Matlab), MSNM 2020, §3.2.2–§3.4.
Recall from Plane Stress Part 1Ingat semula Plane Stress Part 1
Two questions before class. Not graded; one try each.Dua soalan sebelum kelas. Tidak bergred; satu cubaan setiap satu.
The text’s shear strain εxy isShear strain εxy dalam teks ialah
The text’s εxy is the engineering shear strain γxy, the sum of the two slopes; the (1 − ν)/2 entry of [E] is built for it. Row 3 of [B] comes straight from this definition.εxy dalam teks ialah engineering shear strain γxy, iaitu jumlah dua kecerunan; entri (1 − ν)/2 dalam [E] dibina untuknya. Baris 3 [B] terhasil terus daripada definisi ini.
In the 2 × 2 cantilever mesh, global node 5 carries the dofsDalam mesh cantilever 2 × 2, node global 5 membawa dof
Node n owns U2n−1 (x) and U2n (y): node 5 → U9, U10. Taking the node number for the dof number is the commonest assembly slip.Node n memiliki U2n−1 (x) dan U2n (y): node 5 → U9, U10. Mengambil nombor node sebagai nombor dof ialah kesilapan assembly yang paling biasa.
Can-do checklistSenarai semak boleh-buat
Rate yourself after the exercises. “Without notes” is the exam condition and counts in full; “With notes” counts half.Nilai diri anda selepas latihan. “Tanpa nota” ialah keadaan peperiksaan dan dikira penuh; “Dengan nota” dikira separuh.
Write [B] for a rectangle and fill it with numbers at a point; say why each row adds up to zero.Menulis [B] bagi segi empat tepat dan mengisinya dengan nombor pada satu titik; menyatakan mengapa setiap baris berjumlah sifar.
Compute s1…s12 for any a, b, E, ν, t, place them in [k] and check a row.Mengira s1…s12 untuk sebarang a, b, E, ν, t, meletakkannya dalam [k] dan menyemak satu baris.
Turn a uniform or linearly varying edge load and a self-weight into nodal loads in the right dofs.Menukar edge load yang seragam atau berubah secara linear dan self-weight kepada nodal load pada dof yang betul.
Find any global entry KIJ from the dof map, and check a solution by its reactions and by beam theory.Mencari mana-mana entri global KIJ daripada dof map, dan menyemak penyelesaian melalui reaction-nya dan beam theory.
Files and settingsFail dan tetapan
Progress is stored in this browser only. It is your own record, not a gradebook; graded work stays in eLearning.Kemajuan disimpan dalam pelayar ini sahaja. Ia rekod anda sendiri, bukan buku gred; kerja bergred kekal di eLearning.
NotesNota
What to doApa perlu dibuat Read each section, answer its Quick check (not graded, one try), then tick “Mark this section as read” at the end of the section. “Key points only” is the revision view.Baca setiap seksyen, jawab Semak cepat di dalamnya (tidak bergred, satu cubaan), kemudian tanda “Tanda seksyen ini sudah dibaca” di hujung seksyen. “Perkara utama sahaja” ialah paparan ulang kaji.
Counts toward progress:Dikira dalam kemajuan: sections ticked as read.seksyen yang ditanda sudah dibaca.
1. Galerkin and integration by parts, two unknowns per nodeGalerkin dan integration by parts, dua unknown setiap node
Plane Stress Part 1 ended with two results: the PDE [∂][E][∂]ᵀ{u}ᵀ = −{f} and the interpolation {u}ᵀ = [N]{û}ᵀ with the 2 × 8 matrix [N]. Substitute the second into the first and, as in Heat Transfer Part 2, it does not balance: a bilinear displacement field cannot satisfy two second-order equations at every point.Plane Stress Part 1 berakhir dengan dua hasil: PDE [∂][E][∂]ᵀ{u}ᵀ = −{f} dan interpolation {u}ᵀ = [N]{û}ᵀ dengan matrix 2 × 8 [N]. Gantikan yang kedua ke dalam yang pertama dan, seperti dalam Heat Transfer Part 2, ia tidak seimbang: displacement field bilinear tidak mampu memenuhi dua persamaan peringkat kedua pada setiap titik.
Integration by parts, then name the termsIntegration by parts, kemudian namakan term
Exactly as for heat, one [∂] moves from the displacements onto the weight [N]ᵀ (Green’s theorem, applied to each row). The derivative order falls to one, so a bilinear field is admissible, and a boundary integral appears that carries the traction, the natural condition of Plane Stress Part 1.Sama seperti haba, satu [∂] berpindah daripada displacement kepada weight [N]ᵀ (teorem Green, digunakan pada setiap baris). Peringkat derivative turun kepada satu, jadi field bilinear boleh diterima, dan satu integral sempadan muncul yang membawa traction, iaitu natural condition dalam Plane Stress Part 1.
| TermTerm | NameNama | What it holdsApa yang dikandunginya | Heat analogueAnalog haba |
|---|---|---|---|
| [k] | stiffness matrix | material, geometry and thickness of the element; symmetric; 8 × 8bahan, geometri dan ketebalan element; simetri; 8 × 8 | conductance matrix, 4 × 4conductance matrix, 4 × 4 |
| {q} | body-force vectorbody-force vector | forces per unit volume inside the element (self-weight), shared among its nodesdaya per unit isi padu di dalam element (self-weight), dikongsi antara node-nya | source vectorsource vector |
| {b} | traction vectortraction vector | forces on edges with a known traction; zero on free edges; discarded with the equation on fixed edgesdaya pada tepi dengan traction diketahui; sifar pada tepi bebas; dibuang bersama persamaan pada tepi fixed | boundary-flux vectorboundary-flux vector |
Key pointsPerkara utama
- Weight the residual with [N]ᵀ: eight equations per four-node element.Weight residual dengan [N]ᵀ: lapan persamaan setiap four-node element.
- Integration by parts lowers the derivative order and exposes the traction as the natural condition.Integration by parts merendahkan peringkat derivative dan mendedahkan traction sebagai natural condition.
- [k]{û}ᵀ = {q} + {b} = {r}: stiffness, body force, traction; the same structure as heat with two unknowns per node.[k]{û}ᵀ = {q} + {b} = {r}: stiffness, body force, traction; struktur yang sama seperti haba dengan dua unknown setiap node.
Weighting the plane-stress residual with [N]ᵀ gives, for one four-node element,Weight residual plane stress dengan [N]ᵀ memberi, bagi satu four-node element,
[N]ᵀ is 8 × 2 and the residual is 2 × 1, so the weighted residual is an 8 × 1 vector: one equation for each of u1, v1, …, u4, v4. Heat transfer gave four, one per node, because it had one unknown per node.[N]ᵀ ialah 8 × 2 dan residual ialah 2 × 1, jadi weighted residual ialah vector 8 × 1: satu persamaan bagi setiap u1, v1, …, u4, v4. Heat transfer memberi empat, satu setiap node, kerana ia ada satu unknown setiap node.
2. The strain–displacement matrix [B] and the stiffness matrixStrain–displacement matrix [B] dan stiffness matrix
In heat transfer the 2 × 4 matrix {∂}ᵀ{N} turned the nodal temperatures into the gradient. Its plane-stress counterpart turns the eight nodal displacements into the three strains, and it has a name of its own: the strain–displacement matrix [B].Dalam heat transfer, matrix 2 × 4 {∂}ᵀ{N} menukar nodal temperature kepada gradient. Padanannya dalam plane stress menukar lapan nodal displacement kepada tiga strain, dan ia ada nama sendiri: strain–displacement matrix [B].
- Row 1 (εxx = ∂u/∂x): ∂Ni/∂x under each ui, zero under each vi.Baris 1 (εxx = ∂u/∂x): ∂Ni/∂x di bawah setiap ui, sifar di bawah setiap vi.
- Row 2 (εyy = ∂v/∂y): ∂Ni/∂y under each vi, zero under each ui.Baris 2 (εyy = ∂v/∂y): ∂Ni/∂y di bawah setiap vi, sifar di bawah setiap ui.
- Row 3 (εxy = ∂u/∂y + ∂v/∂x, the engineering shear strain): ∂Ni/∂y under ui and ∂Ni/∂x under vi.Baris 3 (εxy = ∂u/∂y + ∂v/∂x, engineering shear strain): ∂Ni/∂y di bawah ui dan ∂Ni/∂x di bawah vi.
Worked step · [B] at (1.5, 0.3) of a 2 m × 1 m element (Activity 1)Langkah bekerja · [B] pada (1.5, 0.3) bagi element 2 m × 1 m (Aktiviti 1)
The derivatives are those of the Heat Transfer Part 1 point: ∂Ni/∂x = {−0.35, 0.35, 0.15, −0.15} and ∂Ni/∂y = {−0.25, −0.75, 0.75, 0.25} m−1. Place them in the three rows:Derivative-nya ialah derivative titik Heat Transfer Part 1: ∂Ni/∂x = {−0.35, 0.35, 0.15, −0.15} dan ∂Ni/∂y = {−0.25, −0.75, 0.75, 0.25} m−1. Letakkannya dalam tiga baris:
Each row adds up to zero, and in row 3 the u-entries and the v-entries add up to zero separately. The reason is physical: a rigid translation (all ui equal, all vi equal) produces no strain. Keep this [B]; Plane Stress Part 3 multiplies it by the nodal displacements to get the stresses.Setiap baris berjumlah sifar, dan dalam baris 3 entri u dan entri v masing-masing berjumlah sifar. Sebabnya fizikal: translation tegar (semua ui sama, semua vi sama) tidak menghasilkan strain. Simpan [B] ini; Plane Stress Part 3 mendarabnya dengan nodal displacement untuk mendapatkan stress.
The element stiffness matrix in [B] formElement stiffness matrix dalam bentuk [B]
Integration by parts left the product [N]ᵀ[∂] on the left and [∂]ᵀ[N] on the right of [E]. Both are [B] (the left one transposed), so the stiffness matrix is the same [B] twice with the material between:Integration by parts meninggalkan hasil darab [N]ᵀ[∂] di kiri dan [∂]ᵀ[N] di kanan [E]. Kedua-duanya ialah [B] (yang kiri di-transpose), jadi stiffness matrix ialah [B] yang sama dua kali dengan bahan di antaranya:
- Thickness. The integrals are per unit thickness; the element has thickness t, so every entry of [k] and every load is multiplied by t.Ketebalan. Integral itu per unit ketebalan; element mempunyai ketebalan t, jadi setiap entri [k] dan setiap load didarab dengan t.
- Symmetry. [B]ᵀ[E][B] is symmetric whenever [E] is, and integration keeps it so: at most 36 distinct entries, and for a rectangle only twelve (§3).Simetri. [B]ᵀ[E][B] simetri apabila [E] simetri, dan integration mengekalkannya: paling banyak 36 entri berbeza, dan bagi segi empat tepat hanya dua belas (§3).
- Meaning. kij is the force at dof i needed to give dof j a unit displacement with every other dof held: the stiffness influence coefficient of structural analysis.Maksud. kij ialah daya pada dof i yang diperlukan untuk memberi dof j satu unit displacement dengan setiap dof lain ditahan: stiffness influence coefficient dalam analisis struktur.
How software integrates [k]Bagaimana perisian meng-integrate [k]

Key pointsPerkara utama
- [B] = [∂]ᵀ[N], 3 × 8: strains from nodal displacements; every row adds up to zero.[B] = [∂]ᵀ[N], 3 × 8: strain daripada nodal displacement; setiap baris berjumlah sifar.
- [k] = t∫[B]ᵀ[E][B] dA, 8 × 8, symmetric; kij = force at i per unit displacement at j.[k] = t∫[B]ᵀ[E][B] dA, 8 × 8, simetri; kij = daya pada i per unit displacement pada j.
- Software integrates at 2 × 2 Gauss points; for a rectangle that is exact.Perisian meng-integrate pada 2 × 2 Gauss point; bagi segi empat tepat ia tepat.
Why must every row of [B] add up to zero?Mengapa setiap baris [B] mesti berjumlah sifar?
Giving every node the same u (or the same v) moves the element without deforming it, so [B]{û}ᵀ must be zero; that is a sum of the u-entries (or the v-entries) of each row. [B] is 3 × 8 and cannot be symmetric; the Kronecker property is true but is not the reason.Memberi setiap node u yang sama (atau v yang sama) menggerakkan element tanpa mengubah bentuknya, jadi [B]{û}ᵀ mesti sifar; itu ialah jumlah entri u (atau entri v) bagi setiap baris. [B] ialah 3 × 8 dan tidak boleh simetri; sifat Kronecker itu benar tetapi bukan sebabnya.
3. The 8 × 8 stiffness matrix and its twelve valuesStiffness matrix 8 × 8 dan dua belas nilainya
One entry in full is enough to see how all of them are made. Take k34: row 3 belongs to u2 and column 4 to v2, so it couples the horizontal and the vertical displacement of the same node.Satu entri secara penuh sudah cukup untuk melihat bagaimana semuanya dibina. Ambil k34: baris 3 milik u2 dan lajur 4 milik v2, jadi ia menggandingkan displacement mendatar dan menegak bagi node yang sama.
k₃₄ in three steps (text (3.39)–(3.41))k₃₄ dalam tiga langkah (teks (3.39)–(3.41))
The integral of −x(b − y)/(a²b²) over the rectangle is −1/4, so k34 does not depend on a or b: the coupling between u and v at a node is the same for every rectangle. The text writes the result as Et/(8(ν − 1)).Integral −x(b − y)/(a²b²) atas segi empat tepat ialah −1/4, jadi k34 tidak bergantung pada a atau b: gandingan antara u dan v pada satu node adalah sama bagi setiap segi empat tepat. Teks menulis hasilnya sebagai Et/(8(ν − 1)).
Twelve values, one patternDua belas nilai, satu pola
Repeating the integration for the other entries gives only twelve distinct values, arranged by the symmetry of the rectangle. Rows 1, 3, 5, 7 belong to u1…u4 and rows 2, 4, 6, 8 to v1…v4; k34 is s6. You do not memorise the pattern: you read it correctly from the text.Mengulang integration untuk entri lain hanya memberi dua belas nilai berbeza, disusun mengikut simetri segi empat tepat. Baris 1, 3, 5, 7 milik u1…u4 dan baris 2, 4, 6, 8 milik v1…v4; k34 ialah s6. Anda tidak perlu menghafal pola ini: anda membacanya dengan betul daripada teks.
Worked example · the cantilever element (Activity 2, text (3.44))Contoh bekerja · element cantilever (Aktiviti 2, teks (3.44))
a = 2 m, b = 1 m, t = 0.01 m, E = 200 × 10⁶ N m−2, ν = 0.33, so 1 − ν² = 0.8911:a = 2 m, b = 1 m, t = 0.01 m, E = 200 × 10⁶ N m−2, ν = 0.33, jadi 1 − ν² = 0.8911:
| 1 · u1 | 2 · v1 | 3 · u2 | 4 · v2 | 5 · u3 | 6 · v3 | 7 · u4 | 8 · v4 | |
|---|---|---|---|---|---|---|---|---|
| 1 · u1 | 0.875 | 0.373 | −0.123 | −0.003 | −0.438 | −0.373 | −0.314 | 0.003 |
| 2 · v1 | 0.373 | 1.622 | 0.003 | 0.623 | −0.373 | −0.811 | −0.003 | −1.434 |
| 3 · u2 | −0.123 | 0.003 | 0.875 | −0.373 | −0.314 | −0.003 | −0.438 | 0.373 |
| 4 · v2 | −0.003 | 0.623 | −0.373 | 1.622 | 0.003 | −1.434 | 0.373 | −0.811 |
| 5 · u3 | −0.438 | −0.373 | −0.314 | 0.003 | 0.875 | 0.373 | −0.123 | −0.003 |
| 6 · v3 | −0.373 | −0.811 | −0.003 | −1.434 | 0.373 | 1.622 | 0.003 | 0.623 |
| 7 · u4 | −0.314 | −0.003 | −0.438 | 0.373 | −0.123 | 0.003 | 0.875 | −0.373 |
| 8 · v4 | 0.003 | −1.434 | 0.373 | −0.811 | −0.003 | 0.623 | −0.373 | 1.622 |
[k] of the cantilever element × 10⁶ N m−1, text (3.44); rows and columns in the dof order u1, v1, …, u4, v4.[k] element cantilever × 10⁶ N m−1, teks (3.44); baris dan lajur mengikut susunan dof u1, v1, …, u4, v4.
All four elements of the 2 × 2 mesh have the same a, b and material, so this one matrix serves the whole mesh. The entries ±0.003 (s4, s8) are small because 1 − 3ν = 0.01 for ν = 0.33.Keempat-empat element mesh 2 × 2 mempunyai a, b dan bahan yang sama, jadi satu matrix ini mencukupi untuk seluruh mesh. Entri ±0.003 (s4, s8) kecil kerana 1 − 3ν = 0.01 bagi ν = 0.33.
Key pointsPerkara utama
- Twelve values s1…s12 in the pattern (3.42); the u–v couplings do not depend on a or b.Dua belas nilai s1…s12 dalam pola (3.42); gandingan u–v tidak bergantung pada a atau b.
- Cantilever element: k11 = s1 = 0.875 × 10⁶, s2 = 0.373 × 10⁶, s9 = 1.622 × 10⁶ N m−1.Element cantilever: k11 = s1 = 0.875 × 10⁶, s2 = 0.373 × 10⁶, s9 = 1.622 × 10⁶ N m−1.
- u-entries and v-entries of every row sum to zero; three zero-energy modes.Entri u dan entri v setiap baris berjumlah sifar; tiga zero-energy mode.
For ν = 1/3, the entries s4 and s8 of a rectangle areBagi ν = 1/3, entri s4 dan s8 bagi segi empat tepat ialah
s8 = −s4 = Et(1 − 3ν)/(8(1 − ν²)), and 1 − 3ν = 0 at ν = 1/3. That is why the cantilever element (ν = 0.33) shows only ±0.003 × 10⁶ there.s8 = −s4 = Et(1 − 3ν)/(8(1 − ν²)), dan 1 − 3ν = 0 pada ν = 1/3. Sebab itulah element cantilever (ν = 0.33) hanya menunjukkan ±0.003 × 10⁶ di situ.
4. Loads: edge tractions and body forcesLoad: traction tepi dan body force
The load vector follows the same rule as in heat transfer, now in the plane-stress dof order: along an edge only the shape functions of its two nodes are non-zero, and along that edge they are the linear functions of a bar. A uniform load is therefore shared equally by the two nodes, in the direction it acts and nowhere else.Load vector mengikut peraturan yang sama seperti heat transfer, kini dalam susunan dof plane stress: di sepanjang satu tepi hanya shape function dua node-nya yang bukan sifar, dan di sepanjang tepi itu ia fungsi linear bar. Load yang seragam oleh itu dikongsi sama rata oleh dua node, dalam arah ia bertindak dan tidak di tempat lain.
A linearly varying edge loadEdge load yang berubah secara linear
If the load per unit length varies from pi at one node to pj at the other, the node at the heavier end takes more. The shares are those of the bar element, and their total and their moment about either node are the two checks:Jika load per unit panjang berubah daripada pi pada satu node kepada pj pada node yang lain, node di hujung yang lebih berat menerima lebih. Bahagiannya sama seperti bar element, dan jumlahnya serta momennya terhadap mana-mana node ialah dua semakan:
Body force: self-weightBody force: self-weight
A uniform body force is shared by all four nodes, a quarter each, because every Ni integrates to ab/4. Self-weight, fy = −ρg = −γ (unit weight γ), acts in −y, so it loads only the v-dofs:Body force yang seragam dikongsi oleh keempat-empat node, satu perempat setiap satu, kerana setiap Ni ber-integral kepada ab/4. Self-weight, fy = −ρg = −γ (berat unit γ), bertindak dalam −y, jadi ia hanya membebankan dof v:
Worked example · the global load vector of the cantilever (Activity 3, text (3.47))Contoh bekerja · global load vector cantilever (Aktiviti 3, teks (3.47))
Elements 2 (nodes 2, 3, 6, 5) and 4 (nodes 5, 6, 9, 8) have their side 2–3 on the free end. Each puts −2500 N on its local v2 and v3, which are global U6, U12 for element 2 and U12, U18 for element 4. Node 6 is shared, so U12 collects −5000 N: the full 5 kN m−1 × 1 m of edge around it.Element 2 (node 2, 3, 6, 5) dan 4 (node 5, 6, 9, 8) mempunyai sisi 2–3 pada hujung bebas. Setiap satu meletakkan −2500 N pada v2 dan v3 local-nya, iaitu U6, U12 global bagi element 2 dan U12, U18 bagi element 4. Node 6 dikongsi, jadi U12 mengumpul −5000 N: keseluruhan 5 kN m−1 × 1 m tepi di sekelilingnya.
The six entries Bi at the fixed dofs 1, 2, 7, 8, 13 and 14 are the unknown support reactions. The total applied load is −10 000 N, which the reactions must balance.Enam entri Bi pada dof fixed 1, 2, 7, 8, 13 dan 14 ialah support reaction yang tidak diketahui. Jumlah load yang dikenakan ialah −10 000 N, yang mesti diimbangi oleh reaction.
Key pointsPerkara utama
- Uniform edge load: p ℓ/2 at each of the two edge nodes, only in the dofs of its direction.Edge load seragam: p ℓ/2 pada setiap satu daripada dua node tepi, hanya pada dof arahnya.
- Linearly varying: ℓ(2pi + pj)/6 and ℓ(pi + 2pj)/6; self-weight: γtab/4 in each v-dof.Berubah secara linear: ℓ(2pi + pj)/6 dan ℓ(pi + 2pj)/6; self-weight: γtab/4 pada setiap dof v.
- Cantilever: R6 = R18 = −2500 N, R12 = −5000 N; six unknown reactions.Cantilever: R6 = R18 = −2500 N, R12 = −5000 N; enam reaction yang tidak diketahui.
The self-weight of an element a × b of thickness t and unit weight γ is shared asSelf-weight bagi element a × b dengan ketebalan t dan berat unit γ dikongsi sebagai
Each shape function integrates to ab/4, so each node carries a quarter of the weight γtab, and the weight acts in −y: only v1…v4 are loaded. Putting it all on the bottom nodes, or into the u-dofs, changes the load the element feels.Setiap shape function ber-integral kepada ab/4, jadi setiap node membawa satu perempat berat γtab, dan berat itu bertindak dalam −y: hanya v1…v4 dibebankan. Meletakkan semuanya pada node bawah, atau pada dof u, mengubah load yang dirasai element.
5. Assembly with two dofs per node, supports and the solutionAssembly dengan dua dof setiap node, support dan penyelesaian
Assembly is the heat-transfer scatter-add with one extra step: node numbers become dof numbers. Node n owns U2n−1 (x) and U2n (y), so each element’s four nodes give eight global dofs, in local order.Assembly ialah scatter-add heat transfer dengan satu langkah tambahan: nombor node menjadi nombor dof. Node n memiliki U2n−1 (x) dan U2n (y), jadi empat node setiap element memberi lapan dof global, mengikut susunan local.
| ElementElement | Nodes (1)(2)(3)(4)Node (1)(2)(3)(4) | Dofs in local order 1…8Dof mengikut susunan local 1…8 |
|---|---|---|
| 1 | 1, 2, 5, 4 | 1 2 · 3 4 · 9 10 · 7 8 |
| 2 | 2, 3, 6, 5 | 3 4 · 5 6 · 11 12 · 9 10 |
| 3 | 4, 5, 8, 7 | 7 8 · 9 10 · 15 16 · 13 14 |
| 4 | 5, 6, 9, 8 | 9 10 · 11 12 · 17 18 · 15 16 |

Worked step · entries of the assembled [K] (Activity 4)Langkah bekerja · entri [K] yang di-assemble (Aktiviti 4)
Local k11 of element 2 (its local node 1 is global node 2) goes to K3,3; local k34 = s6 = −0.373 × 10⁶ goes to K5,6, the u–v coupling of node 3. The diagonals collect s1 once for every element at the node:k11 local bagi element 2 (node local 1-nya ialah node global 2) pergi ke K3,3; k34 local = s6 = −0.373 × 10⁶ pergi ke K5,6, gandingan u–v bagi node 3. Pepenjuru mengumpul s1 sekali bagi setiap element pada node itu:
K9,10 is zero: at an interior node of a regular mesh the u–v couplings of the four elements cancel in pairs. K3,3 = 1.751 and K9,9 = 3.501 (× 10⁶) are visible in the text’s (3.46): a quick check of your dof map.K9,10 sifar: pada node dalaman mesh yang seragam, gandingan u–v empat element saling membatalkan secara berpasangan. K3,3 = 1.751 dan K9,9 = 3.501 (× 10⁶) kelihatan dalam (3.46) teks: semakan pantas untuk dof map anda.
Supports: strike the fixed dofsSupport: buang dof yang fixed
The three nodes of the left edge are fixed in both directions: six known displacements, all zero. As for a prescribed temperature, their rows and columns are deleted; nothing moves to the right-hand side because the known values are zero. Twelve equations remain, text (3.48), which MATLAB solves.Tiga node pada tepi kiri fixed dalam kedua-dua arah: enam displacement yang diketahui, semuanya sifar. Seperti temperature yang ditetapkan, baris dan lajurnya dibuang; tiada yang berpindah ke sebelah kanan kerana nilai yang diketahui ialah sifar. Dua belas persamaan tinggal, teks (3.48), yang diselesaikan oleh MATLAB.
| NodeNode | u (m) | v (m) | WhereDi mana |
|---|---|---|---|
| 1, 4, 7 | 0 | 0 | fixed edgetepi fixed |
| 2 | −0.030966 | −0.045520 | bottom, mid-spanbawah, tengah rentang |
| 5 | 0 | −0.043013 | centrepusat |
| 8 | +0.030966 | −0.045520 | top, mid-spanatas, tengah rentang |
| 3 | −0.041967 | −0.130578 | free end, bottomhujung bebas, bawah |
| 6 | 0 | −0.130421 | free end, mid-heighthujung bebas, tengah ketinggian |
| 9 | +0.041967 | −0.130578 | free end, tophujung bebas, atas |
Read the solution like a structural engineer (Activity 5)Baca penyelesaian seperti jurutera struktur (Aktiviti 5)
- The bottom fibre shortens (u2, u3 < 0) and the top lengthens (u8, u9 > 0) by the same amounts: antisymmetric bending about mid-height, an exact property of this symmetric mesh and load.Gentian bawah memendek (u2, u3 < 0) dan gentian atas memanjang (u8, u9 > 0) dengan jumlah yang sama: lenturan antisimetri tentang tengah ketinggian, sifat tepat mesh dan load yang simetri ini.
- u = 0 at nodes 5 and 6: the mid-height line is the neutral axis.u = 0 pada node 5 dan 6: garis tengah ketinggian ialah neutral axis.
- v is almost the same across the depth of each section (−0.1306, −0.1304, −0.1306 m at the free end): the depth hardly changes, the usual beam assumption.v hampir sama merentasi kedalaman setiap keratan (−0.1306, −0.1304, −0.1306 m di hujung bebas): kedalaman hampir tidak berubah, andaian biasa beam.
Reactions: the first check of any solutionReaction: semakan pertama bagi mana-mana penyelesaian
Evaluate the struck rows, [K]{U} − {R}, at the six fixed dofs. Their sum must balance the applied load and their moment must balance its moment about the fixed edge. The 2 × 2 mesh gives exactly these values, as every correct model must, whatever its mesh:Nilaikan baris yang dibuang, [K]{U} − {R}, pada enam dof fixed. Jumlahnya mesti mengimbangi load yang dikenakan dan momennya mesti mengimbangi momen load itu terhadap tepi fixed. Mesh 2 × 2 memberi nilai ini dengan tepat, seperti yang mesti diberi oleh setiap model yang betul, tanpa mengira mesh-nya:
Key pointsPerkara utama
- Dof map: node n → 2n − 1, 2n; scatter-add 8 × 8 blocks; Kii collects s1 (or s9) once per element at the node.Dof map: node n → 2n − 1, 2n; scatter-add blok 8 × 8; Kii mengumpul s1 (atau s9) sekali bagi setiap element pada node itu.
- Six fixed dofs struck → 12 × 12; tip v = −0.1306 m.Enam dof fixed dibuang → 12 × 12; v hujung = −0.1306 m.
- Reactions: ΣRy = +10 kN, moment 40 kN·m about the fixed edge.Reaction: ΣRy = +10 kN, momen 40 kN·m terhadap tepi fixed.
In the 2 × 2 cantilever, why is u = 0 at nodes 5 and 6?Dalam cantilever 2 × 2, mengapa u = 0 pada node 5 dan 6?
Under bending the fibres above mid-height stretch and those below shorten by the same amount, so the mid-height line neither stretches nor shortens. Only nodes 1, 4 and 7 are supported, and only their rows were struck.Di bawah lenturan, gentian di atas tengah ketinggian meregang dan yang di bawah memendek dengan jumlah yang sama, jadi garis tengah ketinggian tidak meregang mahupun memendek. Hanya node 1, 4 dan 7 yang ditahan, dan hanya barisnya yang dibuang.
6. Is the answer right? Beam theory, refinement and shear lockingAdakah jawapan betul? Beam theory, refinement dan shear locking
The computer will always return numbers. Whether they are right is the engineer’s question, and it has two independent answers: compare with a simpler theory, and refine the mesh until the answer stops moving.Komputer akan sentiasa memulangkan nombor. Sama ada nombor itu betul ialah soalan jurutera, dan ia ada dua jawapan yang bebas: bandingkan dengan teori yang lebih mudah, dan refine mesh sehingga jawapan berhenti berubah.
A hand check: Timoshenko beam theorySemakan tangan: Timoshenko beam theory
The plate is a deep cantilever (L/H = 2) carrying P = qyH = 10 kN at its end. Elementary bending gives PL³/(3EI); because the beam is deep, shear deformation adds PL/(κGA) with κ = 5/6:Plat ini ialah cantilever dalam (L/H = 2) yang membawa P = qyH = 10 kN di hujungnya. Lenturan asas memberi PL³/(3EI); kerana beam ini dalam, shear deformation menambah PL/(κGA) dengan κ = 5/6:

Refine until the answer stops movingRefine sehingga jawapan berhenti berubah
| MeshMesh | dofs | tip v (m)v hujung (m) | vs beam theoryberbanding beam theory |
|---|---|---|---|
| Q4 2 × 2 (handtangan) | 18 | −0.1306 | −32% |
| Q4 4 × 2 | 30 | −0.1643 | −14% |
| Q4 8 × 4 | 90 | −0.1824 | −5.0% |
| Q4 16 × 8 | 306 | −0.1886 | −1.7% |
| Q4 16 × 16 | 578 | −0.1893 | −1.4% |
| Q4 32 × 16 | 1122 | −0.1906 | −0.7% |
| Q4 80 × 80 | 13 122 | −0.1915 | −0.2% |
| Timoshenko beamTimoshenko beam | — | −0.1919 | — |




The same experiment in another code: the error falls with the element size, then levels off near 0.7%, because the load is applied as a uniform end traction while the exact solution assumes a parabolic end shear. Refining removes mesh error, not idealisation error. Images: SolidsPy (J. Gómez and N. Guarín-Zapata), MIT licence.Eksperimen yang sama dalam kod lain: ralat jatuh bersama saiz element, kemudian mendatar dekat 0.7%, kerana load dikenakan sebagai traction hujung yang seragam manakala penyelesaian tepat mengandaikan shear hujung parabolic. Refinement menghapuskan ralat mesh, bukan ralat idealisation. Imej: SolidsPy (J. Gómez dan N. Guarín-Zapata), lesen MIT.
The script, as a diff against Heat Transfer Part 2Skrip, sebagai perbezaan daripada Heat Transfer Part 2
% ---- changes from the heat script ---- Lx = 4; Ly = 2; % plate 4 m x 2 m E = 200e6; nu = 0.33; t = 0.01; Nx = 2; Ny = 2; ndof = 18; % 2 dofs per node a = Lx/Nx; b = Ly/Ny; % 12 closed-form values, text (3.42) s1 = E*t*((1-nu)*a^2 + 2*b^2)/(6*(1-nu^2)*a*b); s2 = E*t*((1-nu) + 2*nu)/(8*(1-nu^2)); % ... s3 to s12 likewise ... k = [s1 s2 s3 s4 s5 s6 s7 s8; ...]; % 8 x 8 pattern of (3.42) % assembly with the dof map (node n -> 2n-1, 2n) K = zeros(ndof); F = zeros(ndof,1); d1 = [1 2 3 4 9 10 7 8]; K(d1,d1) = K(d1,d1) + k; d2 = [3 4 5 6 11 12 9 10]; K(d2,d2) = K(d2,d2) + k; d3 = [7 8 9 10 15 16 13 14]; K(d3,d3) = K(d3,d3) + k; d4 = [9 10 11 12 17 18 15 16]; K(d4,d4) = K(d4,d4) + k; F([6 12 18]) = [-2500 -5000 -2500]; % edge traction, (3.47) fixed = [1 2 7 8 13 14]; free = setdiff(1:ndof, fixed); U = zeros(ndof,1); U(free) = K(free,free) \ F(free); % (3.49)
- Two dofs per node: ndof = 2 × nodes, and each element’s dof list replaces its node list.Dua dof setiap node: ndof = 2 × node, dan senarai dof setiap element menggantikan senarai node-nya.
- The 4 × 4 conductance pattern becomes the 8 × 8 pattern of (3.42).Pola conductance 4 × 4 menjadi pola 8 × 8 daripada (3.42).
- Loads in the v-dofs of the loaded edge; six fixed dofs instead of fixed temperatures. The solve line is unchanged.Load pada dof v tepi yang dibebankan; enam dof fixed menggantikan temperature tetap. Baris penyelesaian tidak berubah.
The Python version, fe_plane_stress.py, is in the Explorer tab; it builds [k] by Gauss quadrature, reproduces (3.49) and runs the convergence study above.Versi Python, fe_plane_stress.py, ada dalam tab Explorer; ia membina [k] melalui Gauss quadrature, menghasilkan semula (3.49) dan menjalankan kajian convergence di atas.
Key pointsPerkara utama
- Two independent checks: beam theory (0.192 m) and mesh refinement; plus the reactions.Dua semakan bebas: beam theory (0.192 m) dan mesh refinement; serta reaction.
- The 2 × 2 Q4 mesh is 32% too stiff: straight-sided elements lock in bending.Mesh Q4 2 × 2 terlalu stiff sebanyak 32%: element bertepi lurus mengalami locking dalam lenturan.
- About 32 elements along the length for 1%; elements long in the bending direction lock most.Kira-kira 32 element sepanjang panjangnya untuk 1%; element yang panjang dalam arah lenturan paling teruk locking.
The same eight Q4 elements as a 4 × 2 mesh or as a 2 × 4 mesh. Which is stiffer?Lapan element Q4 yang sama sebagai mesh 4 × 2 atau sebagai mesh 2 × 4. Yang mana lebih stiff?
2 × 4 gives −0.133 m against −0.164 m for 4 × 2: the 2 m × 0.5 m elements cannot curve along their length, so they develop most spurious shear. The count of elements is the same; their shape is not.2 × 4 memberi −0.133 m berbanding −0.164 m bagi 4 × 2: element 2 m × 0.5 m tidak boleh melengkung sepanjang panjangnya, jadi ia menghasilkan shear palsu paling banyak. Bilangan element sama; bentuknya tidak.
SlidesSlides
What to doApa perlu dibuat Step through with ◀ ▶ (⏮ ⏭ for the first and last), the ← → keys, the mouse wheel or a swipe on the slide, or drag the bar under the slide to any page; the list below jumps to any slide. Download the PDF for printing and offline reading.Gerak dengan ◀ ▶ (⏮ ⏭ untuk slide pertama dan terakhir), kekunci ← →, roda tetikus atau leret pada slide, atau seret bar di bawah slide ke mana-mana halaman; senarai di bawah melompat ke mana-mana slide. Muat turun PDF untuk cetakan dan bacaan luar talian.
Counts toward progress:Dikira dalam kemajuan: each slide you open, once.setiap slide yang anda buka, sekali.
Loading slides…Memuatkan slides…
ExamplesContoh
What to doApa perlu dibuat Model answers in exam form, with suggested marks. Try each question on paper first, then compare step by step; tick it at the end when you have worked through it.Jawapan model dalam bentuk peperiksaan, dengan markah cadangan. Cuba setiap soalan di atas kertas dahulu, kemudian bandingkan langkah demi langkah; tandanya di hujung apabila anda sudah mengikutinya.
Counts toward progress:Dikira dalam kemajuan: examples ticked as worked through.contoh yang ditanda sudah diikuti.
How plane stress is examinedBagaimana plane stress diuji
Model answers in the form the examiner expects: physical picture, plan, numbered steps with every substitution written out, checks, a boxed answer and an interpretation. The marks shown are a suggested scheme; the examiner’s own scheme may differ.Jawapan model dalam bentuk yang pemeriksa jangkakan: gambaran fizikal, pelan, langkah bernombor dengan setiap penggantian ditulis, semakan, jawapan berkotak dan tafsiran. Markah yang ditunjukkan ialah skema cadangan; skema pemeriksa sendiri mungkin berbeza.
| Question typeJenis soalan | Given and wantedDiberi dan dikehendaki | Plan and the usual trapPelan dan perangkap biasa |
|---|---|---|
| Element stiffnessElement stiffness | a, b, E, ν, t; wanted: entries of [k].a, b, E, ν, t; dikehendaki: entri [k]. | s1 to s12, or two integrals. Trap: plane stress and plane strain mixed.s1 hingga s12, atau dua integral. Perangkap: plane stress dan plane strain bercampur. |
| Nodal loadsNodal load | Edge load or body force; wanted: {r}.Edge load atau body force; dikehendaki: {r}. | Equal shares for uniform loads; dofs of the load’s direction only. Trap: traction per area against load per length.Bahagian sama rata bagi load seragam; dof arah load sahaja. Perangkap: traction per luas berbanding load per panjang. |
| Global stiffness entriesEntri global stiffness | Element [k] (often identical elements) and a numbered mesh; wanted: KIJ.[k] element (selalunya element serupa) dan mesh bernombor; dikehendaki: KIJ. | Dof → node → elements containing both → local positions → sum. Trap: node number taken for dof number.Dof → node → element yang mengandungi kedua-duanya → kedudukan local → jumlah. Perangkap: nombor node diambil sebagai nombor dof. |
| Small solutionPenyelesaian kecil | One or two elements, supports, loads; wanted: displacements and reactions.Satu atau dua element, support, load; dikehendaki: displacement dan reaction. | Strike the restrained dofs, use symmetry, solve, check the reactions against the loads.Buang dof yang ditahan, guna simetri, selesaikan, semak reaction berbanding load. |
| ExplainTerangkan | Plane stress against plane strain; why Q4 is stiff in bending; the patch test.Plane stress berbanding plane strain; mengapa Q4 stiff dalam lenturan; patch test. | Claim, reason, consequence — three sentences, in that order.Dakwaan, sebab, akibat — tiga ayat, mengikut urutan itu. |
Notation follows the companion: Ē = E/(1 − ν²), G = Ē(1 − ν)/2, s1…s12 as in text (3.42), nodal loads ri numbered by local dof (u1 = 1, v1 = 2, …). Units kN and m unless stated.Notasi mengikut buku pendamping: Ē = E/(1 − ν²), G = Ē(1 − ν)/2, s1…s12 seperti dalam teks (3.42), nodal load ri dinomborkan mengikut dof local (u1 = 1, v1 = 2, …). Unit kN dan m melainkan dinyatakan.
Example 1 · Stiffness matrix of a square elementContoh 1 · Stiffness matrix bagi element segi empat sama
Physical picture.Gambaran fizikal. k11 is the force in x at node 1 that moves node 1 one unit in x with every other dof held. The movement shortens side 1–2, resisted by Ē = E/(1 − ν²), and shears the element, because node 1 moves while node 4 above it does not, resisted by G. k12, the force in y that comes with it, arises from the Poisson effect and from the same shear.k11 ialah daya dalam x pada node 1 yang menggerakkan node 1 satu unit dalam x dengan setiap dof lain ditahan. Pergerakan itu memendekkan sisi 1–2, ditentang oleh Ē = E/(1 − ν²), dan men-shear element, kerana node 1 bergerak manakala node 4 di atasnya tidak, ditentang oleh G. k12, daya dalam y yang menyertainya, timbul daripada kesan Poisson dan daripada shear yang sama.
Plan. Columns 1 and 2 of [B] → integrate k11 and k12 → s1 to s12 → [k] by the pattern (3.42) → row sums and a rigid rotation.
Step 1 (a) Columns 1 and 2 of [B](a) Lajur 1 dan 2 [B][2]
with Ē = E/(1 − ν²) and G = Ē(1 − ν)/2 in [E].dengan Ē = E/(1 − ν²) dan G = Ē(1 − ν)/2 dalam [E].
Step 2 (a) k₁₁ by integration(a) k₁₁ melalui integration[4]
Step 3 (a) k₁₂ by integration(a) k₁₂ melalui integration[4]
The integral of (a − x)(b − y) over the element is a²b²/4, so k12 contains neither a nor b.Integral (a − x)(b − y) atas element ialah a²b²/4, jadi k12 tidak mengandungi a mahupun b.
Step 4 (b) The twelve values and [k](b) Dua belas nilai dan [k][4]
| 1 · u1 | 2 · v1 | 3 · u2 | 4 · v2 | 5 · u3 | 6 · v3 | 7 · u4 | 8 · v4 | |
|---|---|---|---|---|---|---|---|---|
| 1 · u1 | 5.069 | 1.830 | −3.098 | −0.141 | −2.534 | −1.830 | 0.563 | 0.141 |
| 2 · v1 | 1.830 | 5.069 | 0.141 | 0.563 | −1.830 | −2.534 | −0.141 | −3.098 |
| 3 · u2 | −3.098 | 0.141 | 5.069 | −1.830 | 0.563 | −0.141 | −2.534 | 1.830 |
| 4 · v2 | −0.141 | 0.563 | −1.830 | 5.069 | 0.141 | −3.098 | 1.830 | −2.534 |
| 5 · u3 | −2.534 | −1.830 | 0.563 | 0.141 | 5.069 | 1.830 | −3.098 | −0.141 |
| 6 · v3 | −1.830 | −2.534 | −0.141 | −3.098 | 1.830 | 5.069 | 0.141 | 0.563 |
| 7 · u4 | 0.563 | −0.141 | −2.534 | 1.830 | −3.098 | 0.141 | 5.069 | −1.830 |
| 8 · v4 | 0.141 | −3.098 | 1.830 | −2.534 | −0.141 | 0.563 | −1.830 | 5.069 |
[k] × 10⁶ kN m−1, in the pattern of (3.42); dof order u1, v1, …, u4, v4.[k] × 10⁶ kN m−1, dalam pola (3.42); susunan dof u1, v1, …, u4, v4.
Step 5 (c) Checks(c) Semakan[6]
Check. Symmetric; u- and v-entries of row 1 each sum to zero; a rigid rotation gives zero force. The other rows check in the same way.
Marking scheme (suggested)Skema pemarkahan (cadangan) · 20
| What earns the marksApa yang diberi markah | MarksMarkah |
|---|---|
| Columns of [B] with Ē and G identifiedLajur [B] dengan Ē dan G dikenal pasti | 2 |
| k11: integrand and integrationk11: integrand dan integration | 4 |
| k12: integrand and integrationk12: integrand dan integration | 4 |
| Six distinct values and the square’s equalitiesEnam nilai berbeza dan kesamaan segi empat sama | 4 |
| Row sums and the rigid rotationJumlah baris dan rotation tegar | 6 |
| TotalJumlah | 20 |
Example 2 · Nodal loads, and global entries from a given element matrixContoh 2 · Nodal load, dan entri global daripada matrix element yang diberi
Physical picture.Gambaran fizikal. A uniform load is shared equally by the two nodes of its side, and a uniform body force by the four nodes of the element. A load that grows towards node 3 must put more on node 3, with its resultant at two thirds of the side from node 2. In (d) each global entry collects the stiffness every element contributes between its two dofs; only elements that contain both nodes contribute.Load yang seragam dikongsi sama rata oleh dua node sisinya, dan body force yang seragam oleh empat node element. Load yang membesar ke arah node 3 mesti meletakkan lebih pada node 3, dengan paduannya pada dua pertiga sisi dari node 2. Dalam (d), setiap entri global mengumpul stiffness yang disumbangkan setiap element antara dua dof-nya; hanya element yang mengandungi kedua-dua node menyumbang.
Plan. (a) pℓ/2 → (b) γtab/4 → (c) ℓ(2pi + pj)/6 with total and moment checks → (d) dof → node → elements containing both → local positions → sum.
Step 1 (a) Uniform edge load(a) Edge load seragam[3]
The right side runs from node 2 to node 3, ℓ = 0.5 m; the load is in y, so it enters v2 and v3 (local dofs 4 and 6).Sisi kanan berjalan dari node 2 ke node 3, ℓ = 0.5 m; load dalam y, jadi ia masuk ke v2 dan v3 (dof local 4 dan 6).
Step 2 (b) Self-weight(b) Self-weight[3]
Step 3 (c) Linearly varying load(c) Load yang berubah secara linear[4]
With p2 = 0 and p3 = −25 kN m−1:Dengan p2 = 0 dan p3 = −25 kN m−1:
Step 4 (d) Dofs and connectivity(d) Dof dan connectivity[3]
| 1 · u1 | 2 · v1 | 3 · u2 | 4 · v2 | 5 · u3 | 6 · v3 | 7 · u4 | 8 · v4 | |
|---|---|---|---|---|---|---|---|---|
| 1 · u1 | 2.44 | 0.83 | −1.44 | −0.17 | −1.22 | −0.83 | 0.22 | 0.17 |
| 2 · v1 | 0.83 | 2.44 | 0.17 | 0.22 | −0.83 | −1.22 | −0.17 | −1.44 |
| 3 · u2 | −1.44 | 0.17 | 2.44 | −0.83 | 0.22 | −0.17 | −1.22 | 0.83 |
| 4 · v2 | −0.17 | 0.22 | −0.83 | 2.44 | 0.17 | −1.44 | 0.83 | −1.22 |
| 5 · u3 | −1.22 | −0.83 | 0.22 | 0.17 | 2.44 | 0.83 | −1.44 | −0.17 |
| 6 · v3 | −0.83 | −1.22 | −0.17 | −1.44 | 0.83 | 2.44 | 0.17 | 0.22 |
| 7 · u4 | 0.22 | −0.17 | −1.22 | 0.83 | −1.44 | 0.17 | 2.44 | −0.83 |
| 8 · v4 | 0.17 | −1.44 | 0.83 | −1.22 | −0.17 | 0.22 | −0.83 | 2.44 |
The given element matrix (same local dof order as Example 1), in the units of the question.Matrix element yang diberi (susunan dof local yang sama seperti Contoh 1), dalam unit soalan.
| ElementElement | local 1local 1 | local 2local 2 | local 3local 3 | local 4local 4 |
|---|---|---|---|---|
| 1 | 1 | 2 | 5 | 4 |
| 2 | 3 | 4 | 8 | 7 |
| 3 | 4 | 5 | 9 | 8 |
| 4 | 5 | 6 | 10 | 9 |
[K] is 20 × 20. Global dof 2i − 1 is Ui and 2i is Vi; local dof 2m − 1 is u of local node m and 2m is its v.[K] ialah 20 × 20. Dof global 2i − 1 ialah Ui dan 2i ialah Vi; dof local 2m − 1 ialah u bagi node local m dan 2m ialah v-nya.
Step 5 (d) The four entries(d) Empat entri[7]
- K9,11: U5 and U6. Only element 4 contains nodes 5 and 6, as its local nodes 1 and 2: local dofs 1 and 3.K9,11: U5 dan U6. Hanya element 4 mengandungi node 5 dan 6, sebagai node local 1 dan 2: dof local 1 dan 3.
- K16,16: V8. Node 8 is local node 3 of element 2 (dof 6) and local node 4 of element 3 (dof 8).K16,16: V8. Node 8 ialah node local 3 bagi element 2 (dof 6) dan node local 4 bagi element 3 (dof 8).
- K8,7: V4 and U4. Node 4 is local 4 of element 1, local 2 of element 2 and local 1 of element 3.K8,7: V4 dan U4. Node 4 ialah local 4 bagi element 1, local 2 bagi element 2 dan local 1 bagi element 3.
- K13,20: U7 and V10. Nodes 7 and 10 share no element.K13,20: U7 dan V10. Node 7 dan 10 tidak berkongsi element.
Marking scheme (suggested)Skema pemarkahan (cadangan) · 20
| What earns the marksApa yang diberi markah | MarksMarkah |
|---|---|
| (a) Shares and the dofs that receive them(a) Bahagian dan dof yang menerimanya | 3 |
| (b) Weight and quarter shares(b) Berat dan bahagian satu perempat | 3 |
| (c) Linear shares with total and moment checks(c) Bahagian linear dengan semakan jumlah dan momen | 4 |
| (d) Dof numbering and connectivity table(d) Penomboran dof dan connectivity table | 3 |
| (d) Four entries with the local positions named(d) Empat entri dengan kedudukan local dinamakan | 7 |
| TotalJumlah | 20 |
Example 3 · One element in tension: rollers and a clamped edgeContoh 3 · Satu element dalam tegangan: roller dan tepi clamped
Physical picture.Gambaran fizikal. On rollers the plate is free to contract sideways, so the stress is uniform: it stretches by σ/E = 5 × 10−6 and contracts by νσ/E = 1.5 × 10−6 per metre. That field is linear, and a Q4 contains linear fields, so (a) should come out exact. Clamped, the left side cannot contract, which stiffens the plate near the support: expect a slightly smaller stretch in (b), and the right-hand corners moving towards the centre line.Di atas roller, plat bebas mengecut ke sisi, jadi stress-nya seragam: ia meregang sebanyak σ/E = 5 × 10−6 dan mengecut sebanyak νσ/E = 1.5 × 10−6 per meter. Field itu linear, dan Q4 mengandungi field linear, jadi (a) sepatutnya tepat. Apabila clamped, sisi kiri tidak dapat mengecut, yang menjadikan plat lebih stiff berhampiran support: jangkakan regangan yang sedikit lebih kecil dalam (b), dan bucu kanan bergerak ke arah garis tengah.
Plan. s1 to s12 → loads σtℓ/2 → (a) substitute the exact field into the reduced equations → (b) symmetry u3 = u2, v3 = −v2 → a 2 × 2 system → reactions.
Step 1 Element matrix and loadsMatrix element dan load[4]
For the unit square s9 = s1, s10 = s7, s11 = s5, s12 = s3, s6 = −s2, s8 = −s4.Bagi segi empat sama unit s9 = s1, s10 = s7, s11 = s5, s12 = s3, s6 = −s2, s8 = −s4.
Step 2 (a) Rollers: the exact field satisfies the equations(a) Roller: field tepat memenuhi persamaan[6]
The unknowns are u2, u3, v3, v4. Try the exact field:Unknown-nya ialah u2, u3, v3, v4. Cuba field tepat:
In row u2 (row 3 of [k]) only four terms remain (si in 10⁶ kN m−1, displacements in 10−6 m):Dalam baris u2 (baris 3 [k]) hanya empat term tinggal (si dalam 10⁶ kN m−1, displacement dalam 10−6 m):
Rows u3, v3 and v4 give 5, 0 and 0 in the same way, so the exact field is the solution.Baris u3, v3 dan v4 memberi 5, 0 dan 0 dengan cara yang sama, jadi field tepat itulah penyelesaiannya.
Step 3 (b) Clamped: symmetry gives a 2 × 2 system(b) Clamped: simetri memberi sistem 2 × 2[6]
The unknowns are u2, v2, u3, v3. Plate and load are symmetric about mid-height, so u3 = u2 and v3 = −v2; rows u2 and v2 (rows 3 and 4 of [k]) becomeUnknown-nya ialah u2, v2, u3, v3. Plat dan load simetri tentang tengah ketinggian, jadi u3 = u2 dan v3 = −v2; baris u2 dan v2 (baris 3 dan 4 [k]) menjadi
Step 4 Reactions as a checkReaction sebagai semakan[4]
From rows 1, 2, 7 and 8 of [k]{d} − {r}:Daripada baris 1, 2, 7 dan 8 bagi [k]{d} − {r}:
Check. The horizontal reactions balance the 10 kN pull; the vertical ones are equal and opposite, so they have no resultant.
Marking scheme (suggested)Skema pemarkahan (cadangan) · 20
| What earns the marksApa yang diberi markah | MarksMarkah |
|---|---|
| Element values and the two nodal loadsNilai element dan dua nodal load | 4 |
| (a) Exact field shown to satisfy the reduced equations(a) Field tepat ditunjukkan memenuhi persamaan terkurang | 6 |
| (b) Symmetry, the 2 × 2 system and its solution(b) Simetri, sistem 2 × 2 dan penyelesaiannya | 6 |
| Reactions and the equilibrium checkReaction dan semakan keseimbangan | 4 |
| TotalJumlah | 20 |
ExercisesLatihan
What to doApa perlu dibuat Start the timer. For each question choose Sure or Not sure first, then answer. The summary at the end compares your confidence with your results.Mulakan timer. Untuk setiap soalan pilih Yakin atau Tak pasti dahulu, kemudian jawab. Ringkasan di hujung membandingkan keyakinan anda dengan keputusan.
Counts toward progress:Dikira dalam kemajuan: questions answered (B: checked, or solution viewed).soalan yang dijawab (B: disemak, atau penyelesaian dilihat).
A Concept questionsSoalan konsep · 6 × 1 markmarkah
Click one option. One try: the right answer is marked ✓ and explained.Klik satu pilihan. Satu cubaan: jawapan betul ditanda ✓ dan diterangkan.
The strain–displacement matrix [B] of the four-node plane-stress element isStrain–displacement matrix [B] bagi four-node plane-stress element ialah
[B] = [∂]ᵀ[N] = (3 × 2)(2 × 8) = 3 × 8. The 2 × 8 matrix is [N] itself; the 8 × 8 symmetric matrix is [k] = t∫[B]ᵀ[E][B] dA.[B] = [∂]ᵀ[N] = (3 × 2)(2 × 8) = 3 × 8. Matrix 2 × 8 ialah [N] itu sendiri; matrix simetri 8 × 8 ialah [k] = t∫[B]ᵀ[E][B] dA.
Which entries of a rectangle’s [k] do not depend on the element’s sides a and b?Entri [k] segi empat tepat yang mana tidak bergantung pada sisi element a dan b?
Each u–v coupling integrates the product of an x-derivative and a y-derivative, such as x(b − y)/(a²b²), whose integral over the element is ±1/4 whatever a and b are: s2 = Et/(8(1 − ν)). The diagonals contain a/b and b/a.Setiap gandingan u–v meng-integrate hasil darab satu x-derivative dan satu y-derivative, seperti x(b − y)/(a²b²), yang integral-nya atas element ialah ±1/4 tanpa mengira a dan b: s2 = Et/(8(1 − ν)). Pepenjuru mengandungi a/b dan b/a.
A uniform traction acts in −y on side 2–3 of an element. Its nodal loads go to the local dofsTraction seragam bertindak dalam −y pada sisi 2–3 sesuatu element. Nodal load-nya pergi ke dof local
Along side 2–3 only N2 and N3 are non-zero, and the load acts in y, so it enters the v-dofs of nodes 2 and 3: local dofs 4 and 6, as in text (3.45).Di sepanjang sisi 2–3 hanya N2 dan N3 bukan sifar, dan load bertindak dalam y, jadi ia masuk ke dof v node 2 dan 3: dof local 4 dan 6, seperti dalam teks (3.45).
In the 2 × 2 cantilever mesh, the diagonal entry K9,9 (u of the centre node 5) isDalam mesh cantilever 2 × 2, entri pepenjuru K9,9 (u bagi node pusat 5) ialah
Node 5 belongs to all four elements, and each puts its u–u diagonal s1 = 0.875 × 10⁶ there. K3,3 (node 2, two elements) is 2s1.Node 5 dimiliki oleh keempat-empat element, dan setiap satu meletakkan pepenjuru u–u s1 = 0.875 × 10⁶ di situ. K3,3 (node 2, dua element) ialah 2s1.
The 2 × 2 Q4 mesh gives a tip deflection of 0.131 m; beam theory with shear gives 0.192 m. The coarse mesh isMesh Q4 2 × 2 memberi tip deflection 0.131 m; beam theory dengan shear memberi 0.192 m. Mesh kasar itu
A smaller deflection under the same load means a stiffer model. Bilinear elements cannot curve, so they lock; refining the mesh brings the tip to −0.1906 m on 32 × 16, within 1% of beam theory.Deflection yang lebih kecil di bawah load yang sama bermaksud model yang lebih stiff. Element bilinear tidak boleh melengkung, jadi ia mengalami locking; refinement membawa hujung kepada −0.1906 m pada 32 × 16, dalam 1% daripada beam theory.
After solving the text’s cantilever, the six reactions at the fixed edge mustSelepas menyelesaikan cantilever dalam teks, enam reaction pada tepi fixed mesti
Global equilibrium: the supports carry the whole 10 kN end load and its moment 10 × 4 = 40 kN·m. Zero displacement does not mean zero force; that is exactly where the force is.Keseimbangan global: support membawa keseluruhan load hujung 10 kN dan momennya 10 × 4 = 40 kN·m. Displacement sifar tidak bermaksud daya sifar; di situlah daya itu berada.
B Calculation, new numbers each timePengiraan, nombor baharu setiap kali · 16 marksmarkah
Type your numbers and press Check; each field is marked ✓ or ✗. Retry as often as you like, open hints one at a time, or show the solution. New numbers gives a fresh version.Taip nombor anda dan tekan Semak; setiap ruang ditanda ✓ atau ✗. Cuba semula seberapa kerap, buka petunjuk satu demi satu, atau tunjuk penyelesaian. Nombor baharu memberi versi baharu.
B2 · Nodal loads from a varying edge load and self-weightNodal load daripada edge load yang berubah dan self-weight
4 marksmarkahB3 · Global entries of the 2 × 2 cantilever meshEntri global mesh cantilever 2 × 2
4 marksmarkahB4 · The beam-theory check of a cantilever wallSemakan beam theory bagi dinding cantilever
4 marksmarkahC Find the slipCari kesilapan · 2 marksmarkah
Click the line that contains the error. One try.Klik baris yang mengandungi kesilapan. Satu cubaan.
An element with a = 2 m, b = 1 m and t = 0.01 m carries a downward traction of 500 kN m−2 (force per unit area) on its right side, 2–3. A student finds the nodal loads. One line contains the slip. Which one?Satu element dengan a = 2 m, b = 1 m dan t = 0.01 m membawa traction ke bawah 500 kN m−2 (daya per unit luas) pada sisi kanannya, 2–3. Seorang pelajar mencari nodal load. Satu baris mengandungi kesilapan. Yang mana?
Line 4: the traction acts in y, so the loads belong in the v-dofs of nodes 2 and 3, local dofs 4 and 6: r4 = r6 = −2.5 kN, and r3 = r5 = 0. The total in line 5 could not catch it, because the sum is the same whichever dofs receive the shares; the check that does is to ask which displacement the load does work on.Baris 4: traction bertindak dalam y, jadi load itu milik dof v node 2 dan 3, iaitu dof local 4 dan 6: r4 = r6 = −2.5 kN, dan r3 = r5 = 0. Jumlah pada baris 5 tidak dapat menangkapnya, kerana jumlahnya sama tidak kira dof mana yang menerima bahagian itu; semakan yang berkesan ialah bertanya displacement mana yang dikenakan kerja oleh load itu.
∑ SummaryRingkasan
ExplorerExplorer
What to doApa perlu dibuat Work in the Cantilever Explorer (on a phone, open it in its own tab). Do the 4 tasks on the right; answers to the Explorer’s own questions are under “Check answer” inside it. Python code for assignments and projects, and a panel that runs short scripts in the browser, are at the bottom of this tab.Bekerja dalam Cantilever Explorer (di telefon, buka dalam tab sendiri). Buat 4 tugasan di sebelah kanan; jawapan kepada soalan Explorer ada di bawah “Check answer” di dalamnya. Kod Python untuk assignment dan projek, dan panel untuk menjalankan skrip pendek dalam browser, ada di bahagian bawah tab ini.
Counts toward progress:Dikira dalam kemajuan: tasks done.tugasan yang selesai.
Small screen? Skrin kecil? Open the Explorer in its own tabBuka Explorer dalam tab sendiri. Tasks done there are recorded here when you come back.Tugasan yang selesai di situ direkod di sini apabila anda kembali.
TasksTugasan · Cantilever Explorer
Each task ticks itself when the Explorer sees it done. Tick by hand only if it does not.Setiap tugasan ditanda sendiri apabila Explorer melihatnya selesai. Tanda sendiri hanya jika ia tidak ditanda.
1Press “Text example” (4-node, 2 × 2, the text’s data) and read the tip deflection. Does it match the text’s −0.1306 m, and how far is it from beam theory?Tekan “Text example” (4-node, 2 × 2, data teks) dan baca tip deflection. Adakah ia sepadan dengan −0.1306 m dalam teks, dan sejauh mana ia daripada beam theory?
✓Ticks itself: Ditanda sendiri: detected when the 4-node 2 × 2 mesh is solved with the text’s datadikesan apabila mesh 4-node 2 × 2 diselesaikan dengan data teks ·
2Same eight elements, three layouts: 4-node 4 × 2, 8 × 1 and 2 × 4. Which is stiffest, and why?Lapan element yang sama, tiga susunan: 4-node 4 × 2, 8 × 1 dan 2 × 4. Yang mana paling stiff, dan mengapa?
✓Ticks itself: Ditanda sendiri: detected when the 4-node 4 × 2, 8 × 1 and 2 × 4 meshes have all been solved with the text’s datadikesan apabila mesh 4-node 4 × 2, 8 × 1 dan 2 × 4 semuanya telah diselesaikan dengan data teks ·
3Run the mesh study. Find the coarsest 4-node and 8-node meshes within 2% of beam theory and compare their dofs.Jalankan kajian mesh. Cari mesh 4-node dan 8-node paling kasar dalam 2% daripada beam theory dan bandingkan dof-nya.
✓Ticks itself: Ditanda sendiri: detected when the mesh study runs with the text’s datadikesan apabila kajian mesh dijalankan dengan data teks ·
4Replace the edge traction by a point force at the top corner (same total). What happens at the loaded node, and at mid-span?Gantikan traction tepi dengan point force pada bucu atas (jumlah yang sama). Apa yang berlaku pada node yang dibebankan, dan di tengah rentang?
✓Ticks itself: Ditanda sendiri: detected when the point force at the top corner is solved with the text’s data otherwisedikesan apabila point force pada bucu atas diselesaikan dengan data teks yang lain ·
Python code for assignments and projectsKod Python untuk assignment dan projek
The same steps and notation as the notes, written out in plain Python. Use it to check every matrix and number of a hand calculation, and to run what is impractical by hand: mesh studies, contour plots, your own geometry and loads.Langkah dan notasi yang sama seperti nota, ditulis dalam Python biasa. Guna untuk menyemak setiap matrix dan nombor dalam pengiraan tangan, dan untuk menjalankan apa yang tidak praktikal dibuat dengan tangan: kajian mesh, contour plot, geometri dan beban anda sendiri.
- Unzip, then install Python 3.8 or newer and run
pip install numpy matplotlib(matplotlib only for the plots).Unzip, kemudian pasang Python 3.8 atau lebih baharu dan jalankanpip install numpy matplotlib(matplotlib hanya untuk plot). - Run
python fe_plane_stress.py. It reproduces the course example; README.txt lists the output you should see.Jalankanpython fe_plane_stress.py. Ia menghasilkan semula contoh kursus; README.txt menyenaraikan output yang sepatutnya keluar. - Read the docstring at the top of the file for the notation, numbering and units, and how to call it for your own problem.Baca docstring di bahagian atas fail untuk notasi, penomboran dan unit, dan cara memanggilnya untuk masalah anda sendiri.
A hand calculation is still a hand calculation: in a test or an assignment that asks for the steps, show them yourself.Pengiraan tangan tetap pengiraan tangan: dalam ujian atau assignment yang meminta langkah, tunjukkan sendiri.
Show the codeTunjuk kod · fe_plane_stress.py · 323 linesbaris
"""
SEAA 4243 Finite Element Method - plane stress with 4-node and 8-node elements
==============================================================================
Companion code for the Cantilever and Stress Recovery Explorers (Plane Stress).
It follows the course's steps in the text's notation (Mohd Yassin et al.,
MSNM 2020, Chapter 3):
[E] 3 x 3 material matrix, plane stress or plane strain text (3.10)
[B] strains from nodal displacements, {eps} = [d]^T[N]{u} weak form (3.31)-(3.38)
[k] = t * integral [B]^T [E] [B] dA closed form (3.42);
here by Gauss quadrature (2 x 2 points for the 4-node
element, 3 x 3 for the 8-node), which reproduces (3.42)
loads a uniform traction on the free end becomes consistent nodal
loads: q b/2 to each node (4-node, (3.45)); q b/6, 2q b/3,
q b/6 (8-node, (3.51))
solve assemble (3.46), strike the fixed edge's rows and columns
(3.48), solve (3.49); reactions from the struck rows;
stresses {sigma} = [E][B]{u} at any point; principal stresses,
von Mises, and nodal averaging as in ANSYS.
The built-in problem is the text's cantilever: a plate 4 m x 2 m,
t = 0.01 m, E = 200e6 N/m2, nu = 0.33, fixed at x = 0, loaded at x = 4 by a
uniform downward traction of 5000 N per metre of edge (P = 10 kN in total).
The functions take any L, H, mesh, material, thickness and load.
Mesh and numbering
4-node: nodes left to right, then bottom to top (the text's 2 x 2 mesh:
nodes 1-3 along the bottom). dof 2i-1 = u_i, dof 2i = v_i (1-based).
8-node: corner and mid-side nodes are numbered in the order they are
created; use node_at(res, x, y) to find a node.
Units
Any consistent set. The example uses N and m, so stresses come out in
N/m2 (divide by 1e3 for kN/m2). The traction q is force per unit length
of the loaded edge (N/m), the thickness already included.
Run
python fe_plane_stress.py
reproduces the 2 x 2 cantilever (tip deflection, reactions, stresses
at (3.5, 1.3)), the convergence of the 4- and 8-node meshes against
Timoshenko beam theory, and the growth of the stress at the fixed
corner; saves a von Mises contour plot if matplotlib is installed.
Use for your own model
from fe_plane_stress import solve_cantilever, stress_at, principal, von_mises
res = solve_cantilever(L=4, H=2, Nx=8, Ny=4, E=200e6, nu=0.33, t=0.01, q=-5000, nen=8)
s = stress_at(res, 2.0, 1.3) # (sxx, syy, sxy)
A hand calculation is still a hand calculation: use this script to check your
matrices and numbers, and to run the mesh study a project needs.
Dense numpy solver: fine up to a few thousand degrees of freedom.
Requires Python 3.8+ and numpy; matplotlib only for the plot.
Course SEAA 4243, Faculty of Civil Engineering, Universiti Teknologi Malaysia.
Code: MIT licence.
"""
import numpy as np
# ---------------------------------------------------------------------------
# Material
# ---------------------------------------------------------------------------
def E_matrix(E, nu, plane="stress"):
"""[E] of text (3.10) for plane stress; plane strain when plane = "strain"."""
if plane == "stress":
c = E / (1 - nu ** 2)
return c * np.array([[1, nu, 0], [nu, 1, 0], [0, 0, (1 - nu) / 2]])
c = E / ((1 + nu) * (1 - 2 * nu))
return c * np.array([[1 - nu, nu, 0], [nu, 1 - nu, 0], [0, 0, (1 - 2 * nu) / 2]])
# ---------------------------------------------------------------------------
# Shape functions in natural coordinates xi, eta in [-1, 1]
# ---------------------------------------------------------------------------
XI = np.array([-1, 1, 1, -1, 0, 1, 0, -1], dtype=float)
ETA = np.array([-1, -1, 1, 1, -1, 0, 1, 0], dtype=float)
def shape(nen, xi, eta):
"""N (nen,) and dN/dxi, dN/deta (2, nen) of the 4-node or 8-node element."""
if nen == 4:
N = 0.25 * (1 + XI[:4] * xi) * (1 + ETA[:4] * eta)
dN = np.vstack([0.25 * XI[:4] * (1 + ETA[:4] * eta), 0.25 * ETA[:4] * (1 + XI[:4] * xi)])
return N, dN
N = np.zeros(8); dN = np.zeros((2, 8))
for i in range(4):
cx, cy = XI[i], ETA[i]
N[i] = 0.25 * (1 + cx * xi) * (1 + cy * eta) * (cx * xi + cy * eta - 1)
dN[0, i] = 0.25 * cx * (1 + cy * eta) * (2 * cx * xi + cy * eta)
dN[1, i] = 0.25 * cy * (1 + cx * xi) * (cx * xi + 2 * cy * eta)
N[4] = 0.5 * (1 - xi ** 2) * (1 - eta); dN[:, 4] = [-xi * (1 - eta), -0.5 * (1 - xi ** 2)]
N[5] = 0.5 * (1 + xi) * (1 - eta ** 2); dN[:, 5] = [0.5 * (1 - eta ** 2), -eta * (1 + xi)]
N[6] = 0.5 * (1 - xi ** 2) * (1 + eta); dN[:, 6] = [-xi * (1 + eta), 0.5 * (1 - xi ** 2)]
N[7] = 0.5 * (1 - xi) * (1 - eta ** 2); dN[:, 7] = [-0.5 * (1 - eta ** 2), -eta * (1 - xi)]
return N, dN
def B_matrix(nen, a, b, xi, eta):
"""[B] (3 x 2nen) of an a x b rectangle at (xi, eta): rows exx, eyy, gamma_xy."""
_, dN = shape(nen, xi, eta)
dx, dy = dN[0] * 2 / a, dN[1] * 2 / b
B = np.zeros((3, 2 * nen))
B[0, 0::2] = dx
B[1, 1::2] = dy
B[2, 0::2] = dy
B[2, 1::2] = dx
return B
def element_stiffness(nen, a, b, Em, t):
"""[k] = t * integral [B]^T[E][B] dA by Gauss quadrature (2 x 2 for Q4, 3 x 3 for Q8)."""
g, w = np.polynomial.legendre.leggauss(2 if nen == 4 else 3)
k = np.zeros((2 * nen, 2 * nen))
for xi, wi in zip(g, w):
for eta, wj in zip(g, w):
B = B_matrix(nen, a, b, xi, eta)
k += t * (a / 2) * (b / 2) * wi * wj * B.T @ Em @ B
return k
# ---------------------------------------------------------------------------
# Mesh
# ---------------------------------------------------------------------------
def mesh(nen, L, H, Nx, Ny):
a, b = L / Nx, H / Ny
if nen == 4:
nodes = np.array([(i * a, j * b) for j in range(Ny + 1) for i in range(Nx + 1)], dtype=float)
conn = [[j * (Nx + 1) + i, j * (Nx + 1) + i + 1, (j + 1) * (Nx + 1) + i + 1, (j + 1) * (Nx + 1) + i]
for j in range(Ny) for i in range(Nx)]
return a, b, nodes, np.array(conn, dtype=int)
ids, nodes = {}, []
def node(x, y):
key = (round(x, 9), round(y, 9))
if key not in ids:
ids[key] = len(nodes); nodes.append((x, y))
return ids[key]
conn = []
for j in range(Ny):
for i in range(Nx):
x0, y0 = i * a, j * b
conn.append([node(x0, y0), node(x0 + a, y0), node(x0 + a, y0 + b), node(x0, y0 + b),
node(x0 + a / 2, y0), node(x0 + a, y0 + b / 2), node(x0 + a / 2, y0 + b), node(x0, y0 + b / 2)])
return a, b, np.array(nodes, dtype=float), np.array(conn, dtype=int)
def dofs_of(c):
return np.ravel([[2 * n, 2 * n + 1] for n in c])
# ---------------------------------------------------------------------------
# Solve the cantilever
# ---------------------------------------------------------------------------
def solve_cantilever(L=4.0, H=2.0, Nx=2, Ny=2, E=200e6, nu=0.33, t=0.01, q=-5000.0,
nen=4, plane="stress", load="edge"):
"""Fixed at x = 0; vertical load on x = L. load = "edge": uniform traction q (N/m)
by consistent nodal loads; load = "point": the total q*H at the top-right node."""
Em = E_matrix(E, nu, plane)
a, b, nodes, conn = mesh(nen, L, H, Nx, Ny)
n = len(nodes)
K = np.zeros((2 * n, 2 * n))
F = np.zeros(2 * n)
k = element_stiffness(nen, a, b, Em, t)
for c in conn:
d = dofs_of(c)
K[np.ix_(d, d)] += k # scatter-add by dof map
if load == "point":
F[2 * node_at_xy(nodes, L, H) + 1] += q * H
else:
for c in conn:
if np.isclose(nodes[c[1], 0], L): # element side 2-3 on the loaded edge
if nen == 4:
F[2 * c[1] + 1] += q * b / 2; F[2 * c[2] + 1] += q * b / 2
else:
F[2 * c[1] + 1] += q * b / 6; F[2 * c[5] + 1] += 2 * q * b / 3; F[2 * c[2] + 1] += q * b / 6
fixed_nodes = np.where(np.isclose(nodes[:, 0], 0.0))[0]
known = dofs_of(fixed_nodes)
free = np.setdiff1d(np.arange(2 * n), known)
U = np.zeros(2 * n)
U[free] = np.linalg.solve(K[np.ix_(free, free)], F[free])
Rk = K[known] @ U - F[known] # struck rows: support reactions
Rx, Ry = Rk[0::2].sum(), Rk[1::2].sum()
M = sum(-(nodes[nd, 1] - H / 2) * Rk[2 * i] for i, nd in enumerate(fixed_nodes))
return dict(nen=nen, L=L, H=H, Nx=Nx, Ny=Ny, a=a, b=b, E=E, nu=nu, t=t, q=q, plane=plane, Em=Em,
nodes=nodes, conn=conn, K=K, F=F, U=U, known=known, reactions=dict(Rx=Rx, Ry=Ry, M=M),
tip=U[2 * node_at_xy(nodes, L, H) + 1])
def node_at_xy(nodes, x, y):
return int(np.argmin(np.hypot(nodes[:, 0] - x, nodes[:, 1] - y)))
def node_at(res, x, y):
"""0-based index of the node nearest to (x, y)."""
return node_at_xy(res["nodes"], x, y)
# ---------------------------------------------------------------------------
# Stresses
# ---------------------------------------------------------------------------
def locate(res, x, y):
"""Element containing (x, y) and the natural coordinates of the point in it."""
i = min(res["Nx"] - 1, max(0, int(np.floor(x / res["a"] + 1e-12))))
j = min(res["Ny"] - 1, max(0, int(np.floor(y / res["b"] + 1e-12))))
e = j * res["Nx"] + i
x1, y1 = res["nodes"][res["conn"][e][0]]
return e, 2 * (x - x1) / res["a"] - 1, 2 * (y - y1) / res["b"] - 1
def stress_in_element(res, e, xi, eta):
c = res["conn"][e]
B = B_matrix(res["nen"], res["a"], res["b"], xi, eta)
return res["Em"] @ (B @ res["U"][dofs_of(c)])
def stress_at(res, x, y):
"""(sxx, syy, sxy) at a global point, from the element that contains it."""
e, xi, eta = locate(res, x, y)
return stress_in_element(res, e, xi, eta)
def principal(s):
"""sigma1, sigma2 and the angle (degrees) of sigma1 from x: the Mohr circle."""
c = (s[0] + s[1]) / 2
R = np.hypot((s[0] - s[1]) / 2, s[2])
return c + R, c - R, 0.5 * np.degrees(np.arctan2(2 * s[2], s[0] - s[1]))
def von_mises(s):
"""Plane-stress von Mises stress."""
return np.sqrt(s[0] ** 2 - s[0] * s[1] + s[1] ** 2 + 3 * s[2] ** 2)
def nodal_stresses(res, mode="averaged"):
"""Stresses at the nodes, as post-processors show them.
"unaveraged": each element's own value at its nodes, (ne, nen, 3) - jumps visible;
"averaged": the mean over the elements sharing each node, (n, 3);
"centroid": one value per element, (ne, 3)."""
nen, conn = res["nen"], res["conn"]
if mode == "centroid":
return np.array([stress_in_element(res, e, 0, 0) for e in range(len(conn))])
un = np.array([[stress_in_element(res, e, XI[i], ETA[i]) for i in range(nen)] for e in range(len(conn))])
if mode == "unaveraged":
return un
n = len(res["nodes"])
av = np.zeros((n, 3)); cnt = np.zeros(n)
for e, c in enumerate(conn):
for i, nd in enumerate(c):
av[nd] += un[e, i]; cnt[nd] += 1
return av / cnt[:, None]
# ---------------------------------------------------------------------------
# Beam theory for comparison
# ---------------------------------------------------------------------------
def timoshenko_tip(L, H, E, nu, t, q, plane="stress", kappa=5 / 6):
"""Tip deflection magnitude: bending PL^3/(3EI) plus shear PL/(kappa G A)."""
P = -q * H
I, A, G = t * H ** 3 / 12, t * H, E / (2 * (1 + nu))
Eb = E / (1 - nu ** 2) if plane == "strain" else E
return P * L ** 3 / (3 * Eb * I), P * L / (kappa * G * A)
def plot_von_mises(res, filename="von_mises.png"):
try:
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
except ImportError:
print("matplotlib is not installed: skipping the plot")
return None
xs = np.linspace(0, res["L"], 161); ys = np.linspace(0, res["H"], 81)
X, Y = np.meshgrid(xs, ys)
Z = np.vectorize(lambda x, y: von_mises(stress_at(res, x, y)) / 1e3)(X, Y)
fig, ax = plt.subplots(figsize=(8, 4.6))
# the fixed corners are singular: cap the colour scale at the 99th percentile
levels = MaxNLocator(nbins=20).tick_values(0, np.percentile(Z, 99))
cf = ax.contourf(X, Y, Z, levels=levels, cmap="jet", extend="max")
ax.set_aspect("equal"); ax.set_xlabel("x (m)"); ax.set_ylabel("y (m)")
fig.colorbar(cf, ax=ax, label="von Mises (kN/m2)")
ax.set_title(f"{res['nen']}-node, {res['Nx']} x {res['Ny']} mesh")
fig.savefig(filename, dpi=150, bbox_inches="tight")
plt.close(fig)
print("saved", filename)
return filename
# ---------------------------------------------------------------------------
# Course examples
# ---------------------------------------------------------------------------
if __name__ == "__main__":
np.set_printoptions(precision=4, suppress=True, linewidth=140)
base = dict(L=4.0, H=2.0, E=200e6, nu=0.33, t=0.01, q=-5000.0)
r = solve_cantilever(Nx=2, Ny=2, **base)
print("Text cantilever, 4-node 2 x 2 mesh, (3.44)-(3.49) (nodes 1-3 along the bottom)")
U = r["U"].reshape(-1, 2)
for i, (u, v) in enumerate(U):
print(f" node {i + 1}: u = {u:+.6f} v = {v:+.6f} m")
print(f" tip deflection = {r['tip']:.4f} m (text: -0.131)")
R = r["reactions"]
print(f" reactions: Rx = {R['Rx']:.1f} N, Ry = {R['Ry']:.1f} N, M = {R['M']:.1f} N m (10 kN up, 40 kN m)")
s = stress_at(r, 3.5, 1.3) / 1e3
s1, s2, th = principal(s)
print(f" stress at (3.5, 1.3) = {np.round(s, 1)} kN/m2 (315.2, -44.9, -310.9)")
print(f" principal {s1:.1f}, {s2:.1f} kN/m2 at {th:.1f} deg; von Mises {von_mises(s):.1f} (494.4, -224.1; 636.8)")
bend, shear = timoshenko_tip(**base)
print(f"\nConvergence of the tip deflection; Timoshenko beam {bend:.3f} + {shear:.3f} = {bend + shear:.3f} m")
for nen, meshes in ((4, [(2, 1), (2, 2), (4, 2), (8, 4), (16, 8), (32, 16)]), (8, [(2, 1), (4, 2), (8, 4), (16, 8)])):
for Nx, Ny in meshes:
rr = solve_cantilever(Nx=Nx, Ny=Ny, nen=nen, **base)
print(f" {nen}-node {Nx:2d} x {Ny:<2d} dofs {rr['U'].size:5d} tip {rr['tip']:.4f} m")
r8 = solve_cantilever(Nx=16, Ny=8, nen=8, **base)
print(f"\n8-node 16 x 8: sxx(2.0, 1.3) = {stress_at(r8, 2.0, 1.3)[0] / 1e3:.1f} kN/m2 (beam theory 900);"
f" sxx(3.5, 1.3) = {stress_at(r8, 3.5, 1.3)[0] / 1e3:.1f} (beam 225: end effect, St-Venant)")
print("\nFixed corner: sxx at the centroid of the top-left element grows with every refinement (a singularity)")
for Nx, Ny in ((8, 4), (16, 8), (32, 16)):
rr = solve_cantilever(Nx=Nx, Ny=Ny, **base)
print(f" 4-node {Nx} x {Ny}: {stress_at(rr, rr['a'] / 2, 2 - rr['b'] / 2)[0] / 1e3:.0f} kN/m2")
plot_von_mises(r8, "von_mises_q8_16x8.png")
Run the code in your browserJalankan kod dalam browser trialpercubaan
Edit a short script and run it here, with no installation. Python runs on your own device; nothing is sent anywhere. The first run downloads Python, about 9 MB (about 20 MB once a script plots), and takes 10–30 s; later runs take a second or two.Ubah skrip pendek dan jalankan di sini, tanpa pemasangan. Python berjalan pada peranti anda sendiri; tiada apa yang dihantar ke mana-mana. Larian pertama memuat turun Python, lebih kurang 9 MB (lebih kurang 20 MB apabila skrip membuat plot), dan mengambil 10–30 s; larian seterusnya sesaat dua.
Ctrl+Enter (⌘+Enter) runs. Your edits are kept in this browser; Reset script brings back the original. For an assignment or a project, download the code above and run it on your computer.Ctrl+Enter (⌘+Enter) menjalankan skrip. Suntingan anda disimpan dalam browser ini; Set semula skrip memulihkan yang asal. Untuk assignment atau projek, muat turun kod di atas dan jalankan di komputer anda.
GlossaryGlosari
What to doApa perlu dibuat Look up a term; type in the box to filter. Terms stay in English in both languages; the explanation follows the language you chose.Cari istilah; taip dalam kotak untuk menapis. Istilah kekal dalam bahasa Inggeris dalam kedua-dua bahasa; penjelasan mengikut bahasa yang anda pilih.
Not counted toward progress.Tidak dikira dalam kemajuan.
| Galerkin, vector version | Weighted residual with [N]ᵀ as the weight: eight equations per four-node element, one per dof.Weighted residual dengan [N]ᵀ sebagai weight: lapan persamaan setiap four-node element, satu setiap dof. |
| strain–displacement matrix [B] | [B] = [∂]ᵀ[N], 3 × 8: the strains {ε} = [B]{û} at any point of the element; each row sums to zero.[B] = [∂]ᵀ[N], 3 × 8: strain {ε} = [B]{û} pada mana-mana titik element; setiap baris berjumlah sifar. |
| stiffness matrix [k] | t∫[B]ᵀ[E][B] dA, 8 × 8, symmetric; units force per length (N m−1 or kN m−1).t∫[B]ᵀ[E][B] dA, 8 × 8, simetri; unit daya per panjang (N m−1 atau kN m−1). |
| stiffness influence coefficient kij | Force at dof i needed to give dof j a unit displacement with every other dof held.Daya pada dof i yang diperlukan untuk memberi dof j satu unit displacement dengan setiap dof lain ditahan. |
| s1…s12 | The twelve distinct entries of a rectangle’s [k], text (3.42); s2, s4, s6, s8 (u–v couplings) do not depend on a, b.Dua belas entri berbeza [k] segi empat tepat, teks (3.42); s2, s4, s6, s8 (gandingan u–v) tidak bergantung pada a, b. |
| Ē = E/(1 − ν²) | The plane-stress modulus that multiplies the normal terms of [E]; G = Ē(1 − ν)/2 is the shear modulus.Modulus plane stress yang mendarab term normal dalam [E]; G = Ē(1 − ν)/2 ialah shear modulus. |
| zero-energy (rigid-body) modes | Two translations and one rotation that cause no strain; [k]{d} = 0 for each, so [k] is singular until supported.Dua translation dan satu rotation yang tidak menyebabkan strain; [k]{d} = 0 bagi setiap satu, jadi [k] singular sehingga ditahan. |
| row-sum check | In every row of [k] the u-entries and the v-entries each add up to zero.Dalam setiap baris [k], entri u dan entri v masing-masing berjumlah sifar. |
| Gauss quadrature | Numerical integration at chosen points; 2 × 2 points integrate the Q4 [k] of a rectangle exactly.Numerical integration pada titik terpilih; 2 × 2 titik meng-integrate [k] Q4 segi empat tepat dengan tepat. |
| consistent (equivalent) nodal loads | Nodal forces that do the same work as a distributed load; a uniform edge load gives pℓ/2 to each edge node.Daya nod yang melakukan kerja yang sama seperti load teragih; edge load seragam memberi pℓ/2 kepada setiap node tepi. |
| traction · line load | Traction is force per unit area of edge (kN m−2); a line load is per unit length for the whole thickness (kN m−1) = traction × t.Traction ialah daya per unit luas tepi (kN m−2); line load ialah per unit panjang untuk seluruh ketebalan (kN m−1) = traction × t. |
| body force · self-weight | Force per unit volume; self-weight γ gives −γtab/4 in each v-dof of a Q4 element.Daya per unit isi padu; self-weight γ memberi −γtab/4 pada setiap dof v element Q4. |
| dof map | Node n → dofs 2n − 1 (u) and 2n (v); an element’s eight global dofs in local order.Node n → dof 2n − 1 (u) dan 2n (v); lapan dof global sesuatu element mengikut susunan local. |
| scatter-add | KIJ += kij with I, J from the dof map; the same for {R}.KIJ += kij dengan I, J daripada dof map; sama untuk {R}. |
| bandwidth | How far the non-zeros of [K] spread from the diagonal; nodes numbered close together keep it narrow.Sejauh mana entri bukan sifar [K] tersebar dari pepenjuru; node yang dinomborkan berdekatan memastikannya sempit. |
| fixed dofs · reduced system | Known displacements struck by row–column elimination; 18 − 6 = 12 equations for the 2 × 2 cantilever.Displacement yang diketahui dibuang melalui row–column elimination; 18 − 6 = 12 persamaan bagi cantilever 2 × 2. |
| reactions | The struck rows [K]{U} − {R} at the supports; they must balance the applied loads and their moment.Baris yang dibuang [K]{U} − {R} pada support; ia mesti mengimbangi load yang dikenakan dan momennya. |
| neutral axis | The line that neither stretches nor shortens in bending; mid-height of the cantilever (u = 0 at nodes 5 and 6).Garis yang tidak meregang mahupun memendek dalam lenturan; tengah ketinggian cantilever (u = 0 pada node 5 dan 6). |
| Timoshenko beam | Beam theory with shear deformation: δ = PL³/(3EI) + PL/(κGA), κ = 5/6 for a rectangle; 0.192 m for the text’s cantilever.Beam theory dengan shear deformation: δ = PL³/(3EI) + PL/(κGA), κ = 5/6 bagi segi empat tepat; 0.192 m bagi cantilever dalam teks. |
| shear locking | Bilinear elements with straight edges cannot curve, develop spurious shear in bending and are too stiff; worst when long in the bending direction.Element bilinear bertepi lurus tidak boleh melengkung, menghasilkan shear palsu dalam lenturan dan terlalu stiff; paling teruk apabila panjang dalam arah lenturan. |
| convergence study | The same model on at least three meshes, tabulated, until the result stops changing.Model yang sama pada sekurang-kurangnya tiga mesh, dijadualkan, sehingga hasilnya berhenti berubah. |
| patch test | An element must reproduce a uniform state of stress exactly; one Q4 in uniform tension on rollers does.Element mesti menghasilkan semula keadaan stress seragam dengan tepat; satu Q4 dalam tegangan seragam di atas roller melakukannya. |